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ABSTRACT 

A  general,  modular,  pulse  propagation  model  for  underwater  acoustics  that  is  based  on 
linear  systems  theory  for  sound-speed  profiles  as  a  function  of  depth  is  presented.  The 
development  and  computer  implementation  of  the  model,  together  with  results  from 
preliminary  computer  simulation  studies  involving  the  transmission  of  CW  and  LFM  pulses 
from  a  planar  array  of  complex  weighted  point  sources  is  reported.  The  studies  examined 
free-space  propagation  problems  (i.e.,  no  boundaries)  in  homogeneous  and 
inhomogeneous  media  using  a  transfer  function  of  the  ocean  medium  based  on  the  WKB 
approximation.  The  two  main  outputs  from  the  model  are  the  predicted  complex  acoustic 
field  as  a  function  of  frequency  and  spatial  location  and  the  time-domain  output  electrical 
signal  from  each  element  in  a  receive  planar  array. 
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I.   INTRODUCTION 

For  small-amplitude  acoustic  signals  it  can  be  shown  that  the  linear  acoustic  wave 
equation  governs  the  propagation  of  sound  in  a  fluid  medium  [Ref.  l:pp.  23-31]  [Ref. 
2:pp.  98-105].  Due  to  the  linear  nature  of  this  equation,  the  ocean  medium  can  be  treated 
in  general,  as  a  linear,  time-variant,  space-variant,  random  filter  or  communication  channel 
completely  characterized  by  a  transfer  function  [Ref.  3:pp.  1-6].  As  a  result,  linear 
systems  theory  can  be  applied  to  the  ocean  medium  propagation  problem  and  through  the 
use  of  coupling  equations  [Ref.  4]  analytical  expressions  can  be  derived  for  the  complex 
acoustic  field  and  the  output  electrical  signal  at  each  element  in  a  hydrophone  array  in  terms 
of: 

•  the  frequency  spectrum  of  the  transmitted  pulse 

the  far-field  directivity  function  or  beam  pattern  of  the  transmit  array 

•  the  ocean  medium  transfer  function 

•  the  far-field  directivity  function  of  the  receive  array. 

These  expressions  are  developed  for  planar  arrays,  since  they  include  linear  arrays  and 
single  transducers  as  special  cases,  and  are  amenable  to  computer  simulation  studies. 

The  application  of  Fourier  analysis  and  linear  systems  concepts  to  problems  in  physics 
is  not  unique  and  has  been  applied  extensively  to  the  study  of  optics  [Ref.  5].  What  is 
unique  is  the  application  of  these  ideas  to  the  physics  of  ocean  acoustics  and  the  treatment 
of  the  ocean  medium,  in  general,  as  a  time-variant,  space-variant  linear  system  [Ref.  3:pp. 
7-28]  [Ref.  6]. 


In  this  thesis,  we  will  present  a  general,  modular,  pulse  propagation  model  for 
underwater  acoustics  that  is  based  on  linear  systems  theory  and  the  coupling  equations  for 
sound-speed  profiles  that  are  a  function  of  depth.  Since  the  coupling  equations  are  the 
formal  solution  of  the  pulse  propagation  problem,  and  since  the  coupling  equations  depend 
on  the  transfer  function  of  the  ocean  medium,  the  only  thing  that  changes  from  problem  to 
problem  from  an  ocean  acoustics  point  of  view,  is  the  ocean  medium  transfer  function. 
Therefore,  regardless  of  the  problem  under  consideration,  the  coupling  equations  only  need 
to  be  programmed  once.  Preliminary  computer  simulation  studies  can  then  carried  out  for 
the  specific  case  of  a  free-space  propagation  problem  using  a  transfer  function  based  on  the 
WKB  approximation. 

The  evaluation  of  the  ocean  medium  transfer  function  is  based  on  the  full-wave  solution 
technique,  which  involves  the  evaluation  of  wave  propagation  vector  component  integrals. 
This  method  is  different  from,  but  can  be  related  to,  ray  acoustics  and  normal  mode 
approaches  to  solving  the  propagation  problem.  As  a  consequence,  our  model  of  the  pulse 
propagation  problem  is  termed  a  full-wave  solution.  The  computer  implementation  of  a 
full-wave  solution  to  the  propagation  problem  is  not  new  [Ref.  7]  but  an  implementation  of 
a  model  also  based  on  linear  systems  concepts  is. 

For  the  purposes  of  this  thesis  the  transmit  and  receive  arrays  were  assumed  to  be 
motionless  which  meant  that  we  were  able  to  treat  the  ocean  medium  as  a  linear,  time- 
invariant,  space-variant  filter.  Chapter  2  introduces  and  applies  the  linear  systems 
approach  to  the  time-invariant  ocean  medium  problem  and  develops  the  coupling  equations 
as  far  as  possible  without  making  any  further  simplifying  assumptions  about  the  ocean 
medium.  As  a  consequence  the  results  obtained  can  be  applied  to  the  propagation  of  an 
acoustic  signal  in  any  fluid  medium.  Chapter  3  then  develops  the  coupling  equations  for 
the  case  of  an  ocean  medium  that  is  axially  symmetric  about  the  depth  axis.     This 


development  was  based  on  the  fact  that  axial  symmetry  was  due  to  the  speed  of  sound 
being  a  function  of  depth  only,  i.e.,  the  case  we  wish  to  solve  for.  The  results  of  this 
analysis  can  then  be  applied  to  any  ocean  medium  description  that  exhibits  axial  symmetry. 
Chapter  4  then  goes  on  to  develop  the  ocean  medium  transfer  function  for  the  specific  case 
of  the  ocean  modelled  with  the  WKB  approximation. 

At  this  point,  the  development  of  the  pulse  propagation  model  is  complete  and  is 
capable  of  producing  two  major  outputs.  The  first  output  provides  the  magnitude  and 
phase  of  the  complex  acoustic  field  incident  upon  the  receive  array  as  a  function  of 
frequency  and  spatial  coordinates.  This  information  is  important,  for  example,  for  target 
localization  using  matched-field  techniques.  The  second  output  is  the  time-domain  output 
electrical  signal  at  each  element  in  the  receive  array.  This  information  is  important,  for 
example,  to  illustrate  pulse  distortion  due  to  the  effects  of  dispersion  in  a  waveguide,  and 
as  input  to  space-time  signal  processing  algorithms. 

An  integral  part  of  the  model  was  the  development,  in  Chapter  5,  of  a  signal  generator 
to  simulate  the  transmitted  electrical  pulses  that  would  be  converted  to  acoustic  energy  for 
the  propagation  problem.  It  was  assumed  that  the  pulses  were  narrowband  amplitude-and 
angle-modulated  carriers.  The  signal  generation  scheme  was  based  on  a  truncated  Fourier 
series  and  complex  envelope  representation  of  narrowband  signals. 

Chapters  6  and  7  contain  computer  simulation  results  for  the  specific  case  of  the 
ocean  medium  that  is  characterized  by  the  WKB  approximation.  These  chapters  deal  with 
simulation  results  for  an  ocean  medium  characterized  by  a  constant  speed  of  sound  and  a 
sound-speed  profile  with  a  single  constant  gradient,  respectively. 


II.   LINEAR  SYSTEMS 

A.   THE  LINEAR  SYSTEMS  APPROACH 

Treating  acoustic  signals  as  small  fluctuating  disturbances  in  a  fluid  medium,  it  can  be 
shown  that  the  principles  of  fluid  mechanics  can  be  used  to  derive  an  equation  governing 
the  propagation  of  sound  in  a  fluid  medium.   The  linear  acoustic  wave  equation,  given  by 


2  i       3Z 

V  <p(t,r)-  — -9(t,r)  =  xM(t,r),  (2.1) 

c  (r)  at2 


(where  cp(t,r)  is  the  velocity  potential,  in  m2  /  sec,  at  time  t  and  spatial  location  r,  xM  (t,r) 
is  the  input  acoustic  signal  to  the  medium,  in  1  /  sec,  and  c(r)  is  the  speed  of  sound  in  the 
medium,  in  m  /  sec)  can  be  derived  by  making  a  number  of  restrictive  assumptions  and 
applying  them  to  the  equations  of  state,  continuity  and  momentum  for  a  fluid.  These 
restrictions  allow  us  to  derive  a  simple  equation,  of  the  form  given  by  Eq.  (2.1),  which 
nonetheless  is  adequate  to  describe  most  commonly  encountered  acoustic  phenomena. 
[Ref.  l:pp.  23-31]   [Ref.  2:pp.  98-105]  [Ref.  3:pp.  1-3] 

The  more  commonly  encountered  acoustic  quantities  of  acoustic  pressure,  p(t,r),  in 
N  /  m2  and  the  acoustic  fluid  particle  velocity,  u(t,r),  in  m  /  sec  can  be  obtained  from  the 
velocity  potential  as  follows  : 


p(t,r)  =  -p0^q>(t,r)  (2.2) 


where  p0  is  the  constant  equilibrium  density  of  the  fluid  in  kg  /  m3,  and 


u(t,r)  =  Vcp(t,r).  (2.3) 

It  should  be  noted  here  that  the  concept  of  a  scalar  velocity  potential  arises  from  one  of  the 
restrictions  used  in  the  development  of  Eq.  (2.1),  namely  that  the  particle  motions 
associated  with  sound  waves  are  irrotational.  Therefore,  the  relationships  given  by  Eqs. 
(2.2)  and  (2.3)  are  governed  by  this  assumption. 

The  linear  acoustic  wave  equation  given  by  Eq.  (2.1)  is  so  called  because  it  is  a  linear 
partial  differential  equation  relating  the  velocity  potential  and  input  signal.  Since  Eq.  (2.1) 
can  be  thought  of  as  a  linear  operation  denoted  by 


>  =  v2-^,  (2.4) 

c  (r)  at2 


that  relates  the  velocity  potential  and  input  signal,  we  can  use  a  linear  systems  theory 
approach  to  describe  the  effect  of  the  fluid  medium  on  the  propagation  of  an  input  signal 
through  the  medium. 

Using  a  linear  systems  approach,  we  describe  the  fluid  medium  as  being  a  linear  system 
that  is  characterized  by  an  impulse  response  that  is  a  function  of  both  time  and  space.  To 
then  find  the  output  (i.e.,  the  acoustic  velocity  potential  cp(t,r))  from  such  a  system  due  to 
some  input  (i.e.,  the  acoustic  signal  xM(t,r)),  we  simply  convolve  the  input  with  the 
impulse  response  of  the  system.  The  basic  geometry  of  this  problem  is  illustrated  in 
Fig.  2.1. 

This  operation  can  also  be  done  in  the  frequency  ( in  Hz  )  and  spatial  frequency  (in 
cycles  /  m)  domains  using  the  transfer  function  of  the  system.    Since  this  is  analogous  to 


Fig.  2.1   Ocean  Medium  Geometry 


the  concept  of  filtering,  we  say  that  our  system  is  being  modelled  as  a  time-variant,  space- 
variant,  linear  filter.  [Ref.  3:pp.  3-6] 

Up  until  this  point,  we  have  been  talking  about  any  general  fluid  medium.  We  now 
wish  to  consider  this  problem  with  specific  reference  to  the  ocean  medium.  For  the 
purposes  of  this  thesis,  we  will  model  the  ocean  medium  as  a  time-invariant,  space-variant 
linear  filter.  Time  invariance  means  that  we  will  not  be  considering  Doppler  shifts  between 
the  input  and  output  of  the  filter.  In  relation  to  the  ocean  medium  we  are  modelling,  this 
means  that  we  are  ignoring  such  things  as  motion  between  the  transmitter  and  receiver, 
target  motion  and  discrete  point  scatterers. 

Space  variance  means  that  there  is  a  shift  in  spatial  frequency  between  the  transmitted 
and  received  signal.  This  means  that  sound  speed  is  a  function  of  position  and  there  exists 
a  non-uniform  sound-speed  profile.     As  a  result,  a  transmit  signal  will  scatter  due  to 


refraction.  If  we  think  of  a  ray  interpretation  of  acoustic  propagation,  the  rays  will  bend 
rather  than  travel  in  straight  lines  as  is  the  case  for  homogeneous  media. 

B.    THE  LINEAR,  TIME-INVARIANT,  SPACE- VARIANT  FILTER 

Consider  a  system  modelled  by  a  linear,  time-invariant,  space-variant  filter.  The  time- 
invariant,  space-variant  impulse  response  of  the  filter,  h(x,rQ;r),  will  then  completely 
characterize  the  system  being  modeled.  The  relationship  between  the  input  x(t,r)  and  the 
output  y(t,r)  of  such  a  filter  is  given  by 

y(t,r)=  J         f        x(t-x,r-r0)h(x,r0;r)dxdr0  (2.5) 

where  the  function  h(x,r0;r)  is  the  response  of  the  filter  positioned  at  spatial  location 
r  =  (x,y,z),  due  to  an  impulse  applied  at  position  r0  =  (x0  ,y0  ,z0 ),  x  seconds  ago.  For 
the  specific  case  of  modeling  the  ocean  medium,  the  terms  y(t,r)  and  x(t,r)  in  Eq.  (2.5) 
would  equate  to  cp(t,r)  and  xM  (t,r)  in  Fig.  2.1,  respectively. 

As  an  input  to  our  filter,  consider  a  signal  x(t,r)  of  the  following  form: 


x(t,r)  =  ejMle-j2™-r  (2.6) 


where  f  is  the  input  frequency  in  Hz  and  v  =  (fx.fy^z)  are  tne  inPut  sPatial  frequency 
components  in  cycles  /  meter.  The  spatial  frequencies  are  a  function  of  direction  and 
wavelength  given  by 

fx  =uA  ,  (2.7a) 


=  vA  (2.7b) 


and 


fz=wA  (2.7c) 

where  u,  v  and  w  are  the  direction  cosines  with  respect  to  the  positive  X,  Y  and  Z  axes. 
The  surfaces  of  constant  phase,  commonly  referred  to  as  wavefronts,  for  the  signal 
described  by  Eq.  (2.6)  are  given  by  27tvr  =  constant,  i.e.  a  plane.,  Consequently,  Eq. 
(2.6)  describes  a  time-harmonic  plane  wave  traveling  in  the  direction  of  the  propagation 
vector  k  =  2rcv,  normal  to  the  wavefront.[Ref.  2:  pp.  107-108] 

To  find  the  filter  output  due  to  a  time-harmonic  plane  wave  input,  we  combine 
Eqs.  (2.5)  and  (2.6)  which  yields 


y(t,r)  =  ej27rfte-j27CV,rH(f,v;r),  (2.8) 


where 


H(f,v;r)=f        (      h(T,r0;r)e"j2nfV27W'r°dTdr0  (2.9) 


is  the  transfer  function  of  the  linear,  time-invariant,  space-variant  system  (or  filter).  It  is 
completely  analogous  to  the  transfer  function  H(f)  of  linear  time-invariant  (LTI)  systems 
often  encountered  in  electrical  circuit  theory  and  signal  analysis. 


Equation  (2.9)  is  in  the  form  of  a  multidimensional  Fourier  transform  of  the  impulse 
response  and  can  be  written  as 

H(f,v;r)  =  FTFro{h(T,r0;r)}  (2.10) 

where  FT  represents  the  forward  time-domain  Fourier  transform  given  by 

H(f,r0;r)  =  FT  {  h(x,r0;r) }  =  f      h(x,r0;r)  e'j27cfx  dx  ,  (2.11) 

and  Fr   represents  the  forward  spatial  Fourier  transform  given  by 


H(x,v;r)  =  Fro{  h(x,r0;r)  }  =  j 


h(x,r0;r)  ej27CV'r°  dr0  .  (2.12) 


Comparing  Eqs.  (2.8)  and  (2.9),  we  see  that  we  have  two  methods  of  obtaining  the 
transfer  function  of  our  filter.  Since  the  transfer  function  completely  characterizes  the 
behavior  of  the  filter  and,  hence,  the  ocean  medium  that  the  filter  is  modeling,  we  need 
some  method  of  determining  this  quantity.  Equation  (2.9)  is  not  much  use  for  our 
problem  since  it  assumes  that  we  already  know  the  impulse  response  of  the  filter  which  we 
do  not.  However,  Eq.  (2.8)  is  much  more  useful.  It  says  that  if  we  apply  a  plane  wave 
input  to  the  system,  then  there  is  a  simple  relationship  between  the  input  and  output  that 
determines  the  value  of  the  transfer  function.  However,  this  relationship  is  more  versatile 
than  may  first  appear.  It  tells  us  that  if  we  can  express  an  output  due  to  any  type  of  input 
as  the  product  of  the  time-harmonic  plane  wave  given  by  Eq.  (2.6)  and  some  other  factor, 
then  this  other  factor  is  the  transfer  function. 


Expressing  the  impulse  response  term  in  Eq.  (2.5)  in  terms  of  the  two  dimensional 
inverse  Fourier  transform  of  its  transfer  function  and  combining  with  the  input  term,  it  can 
be  shown  that  the  two  dimensional  Fourier  transform  of  the  output  of  a  linear  time- 
invariant,  space-variant  filter,  given  by  Eq.  (2.5),  has  the  form 

Y(f,P)=  f       J       X(f,v)H(f,v;r)  e"j2*v"P)'r°  dv  dr  (2.13) 


which  is  the  frequency  and  angular  spectrum  of  the  output.  Note  that  the  output  spatial 
frequency  term  is  different,  in  general,  from  that  at  the  input,  and  that  Y  *  XH  .  The 
exception  to  this  rule  is  when  the  ocean  medium  is  homogeneous.  In  this  case  the  filter 
model  is  space-invariant  and  Eq.  (2.13)  collapses  to  the  form  Y  =  XH  ,  which  is  the  same 
as  the  familiar  result  encountered  for  one-dimensional  LTI  systems. 

We  can  determine  the  impulse  response  once  we  know  the  transfer  fountain  by  using 
the  inverse  transform  given  by 

h(x,r0;r)=  f°°    f°°  H(f,v;r)  ej2*Vj27CV'r°  df  dv  ,  (2.14) 

which  can  be  expressed  in  transform  notation  as 


h(x,r0;r)  =  Ff1F;1  {  H(f,v;r)  }  (2  15) 


where  Ff  represents  the  inverse  time-domain  Fourier  transform  and  F~v  represents  the 
inverse  spatial  Fourier  transform.  The  form  of  Ff  and  Fv  is  the  same  as  for  Eqs.  (2. 1 1) 
and  (2.12)  with  the  sign  of  the  exponent  changed  and  the  variable  of  integration  changed  to 
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f  and  v,  respectively.  This  is  the  same  as  for  the  familiar  one-dimensioned  time-frequency 
Fourier  transform. 

Note  that  the  development  in  this  section  has  so  far  been  completely  general  and  can  be 
applied  to  any  linear,  time-invariant,  space-variant  system.  We  have  not  made  any 
assumptions  regarding  the  space  variance  of  the  system,  i.e.,  we  have  not  specified  a 
particular  form  for  the  sound  speed-profile  as  a  function  of  three-dimensional  space. 

C.   THE  PROBLEM  STATEMENT 

Let  us  now  consider  the  underwater  acoustic  problem  from  the  viewpoint  of  linear 
systems  theory.  Our  aim  is  to  be  able  to  compute  the  output  electrical  signal  from  a  receive 
transducer  due  to  an  applied  electrical  signal  at  a  transmit  transducer.  Simply  stated,  this 
process  consists  of  three  steps  : 

•  conversion  of  an  electrical  signal  into  a  transmitted  acoustic  signal, 

•  propagation  of  an  acoustic  signal  through  the  ocean  medium,  and 
conversion  of  a  received  acoustic  signal  into  an  electrical  signal 

Our  intention  is  to  model  each  of  these  steps  as  a  time-invariant,  space-variant  /  invariant 
filter  as  appropriate.  The  system  model  conforming  to  this  plan  is  shown  in  Fig.  2.2  in 
block  diagram  form. 

The  remainder  of  this  chapter  will  deal  with  the  development  of  expressions  describing 
the  behavior  of  the  transmit  and  receive  arrays,  and  the  relationship  or  coupling  equation 
that  relates  the  output,  y(t,r),  to  the  input,  x(t,r).  The  development  of  the  transfer  function 
that  characterizes  the  ocean  medium  will  be  discussed  in  the  next  chapter. 


x(t,r) 

AT(f,r) 
DT(f,v) 

xM(t,r) 

hM(x,r0;r) 
HM(f,v;r) 

yM(t>r) 

AR(f,r) 
DR(f,(3) 

y(t,r) 

X(f,a) 

XM(f>v) 

YM(f,P) 

Y(f,7) 

transmit 
array 

ocean 
medium 

receive 
array 

Fig.  2.2  System  Model 

D.   COUPLING  THE  TRANSMITTED  SIGNAL  TO  THE  MEDIUM 

Consider  now  the  specific  problem  of  transmitting  an  electrical  signal  through  the  ocean 
medium.  The  electrical  signal  will  be  used  to  drive  some  transducer  (or  array  of 
transducers)  which  will  convert  the  electrical  energy  into  acoustic  energy.  We  therefore 
need  to  find  a  relationship  that  describes  this  transformation  of  an  electrical  signal  x(t,r) 
into  an  acoustic  signal  xM  (t,r) .  The  acoustic  signal  represents  the  rate  at  which  fluid 
volume  is  added  at  time  t  and  position  r  per  unit  volume  of  fluid. 

Suppose  we  model  an  infinitesimal  region  of  some  arbitrarily  shaped  transducer  as  a 
linear  filter  with  a  corresponding  impulse  response.  The  output  from  such  a  transducer  at 
a  given  position  r  would  be  given  by  the  time-domain  convolution  of  the  corresponding 
impulse  response  with  the  input  signal.  Translated  into  the  frequency  domain  via  a  Fourier 
transform  with  respect  to  time,  this  relationship  can  be  expressed  as 


XM(f,r)  =  X(f,r)AT(f,r) 


(2.16) 


where  AT  (f,r)  is  the  complex  frequency  response  (or  aperture  function)  of  the  transducer 
at  spatial  location  r.  Taking  the  spatial  Fourier  transform  of  both  sides  of  Eq.  (2. 16),  it 
can  be  shown  that : 
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t(f,v)=  J 


X(f,a)DT(f,v-a)da  (2.17) 


where  DT  (f,a)  is  the  spatial  Fourier  transform  of  the  complex  frequency  response  of  the 
transducer  and  is  known  as  the  far-field  directivity  function  or  beam  pattern.  Note  that  the 
result  given  by  Eq.  (2.17)  is  valid  for  a  completely  arbitrary  volume  aperture. 

If  we  make  the  simplifying  assumption  that  an  identical  electrical  signal  is  applied  at  all 
spatial  locations,  r,  of  the  transducer  then 

x(t,r)  =  x(t).  (2.18) 

Taking  both  the  time  and  spatial  domain  Fourier  transforms  of  a  signal  of  the  form  of 
Eq.(2.18)  and  substituting  into  Eq.  (2.17)  yields 

XM(f,v)  =  X(£)D1(f,v).  (2.19) 

If  the  transmit  aperture  is  a  planar  array,  then  for  most  practical  situations  we  can  consider 
that  an  identical  input  electrical  signal  is  applied  to  all  the  elements  of  the  array  before  any 
electronics  used  for  beamsteering.  As  a  result,  the  restriction  imposed  by  Eq.  (2.18)  does 
not  impose  any  limitations.  What  we  gain,  however,  when  comparing  Eqs.  (2.17)  and 
(2.19),  is  a  significant  simplification  in  our  analysis.  [Ref.  4:p.  4] 

E.    THE  TRANSMIT  ARRAY 

Assume  that  the  transmit  aperture  is  a  planar  array  of  identical  transducer  elements. 
The  product  theorem  for  planar  arrays,  given  by 
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D(f,v)  =  S(v)E(f,v),  (2.20) 

states  that  the  far-field  directivity  function  of  such  an  array  is  equal  to  the  far-field 
directivity  function  of  one  of  the  transducer  elements,  E(f,v) ,  multiplied  by  the  far-field 
directivity  function  of  an  equivalent  array  of  complex  weighted  point  sources,  S(v) .  We 
will  assume  for  the  purposes  of  this  thesis  that  our  transmit  array  is  located  in  the  XY  plane 
and  consists  of  MT  x  NT  (odd)  complex-weighted  point  sources.  We  will  further  assume 
that  the  array  elements  are  equally  spaced  in  the  X  and  Y  directions;  however,  it  should  be 
noted  that  this  is  not  a  constraint  imposed  by  the  use  of  the  product  theorem. 

Since  the  spatial  Fourier  transforms  deal  with  position  relative  to  some  source,  it 
makes  sense  to  locate  the  transmit  array  in  a  convenient  position  to  simplify  the  analysis. 
Accordingly,  we  position  the  array  such  that  its  center  forms  the  origin  of  the  coordinate 
frame  of  reference  in  the  X  and  Z  directions  while  the  ocean  surface  forms  the  center  in  the 
Y  direction,  i.e.,  the  array  is  centered  at  (x0  =  0,yo  =  y-r»zo  =  0)  •  The  geometrical 
interpretation  of  this  configuration  is  shown  in  Fig.  2.3,  which  is  also  the  physical 
representation  of  the  block  diagram  of  Fig.  2.2. 

Given  these  conditions,  Ziomek  [Ref.  3:p.  126]  shows  that  the  far-field  directivity 
function  of  the  transmit  array  is  given  by 


MT'  NT 

DT(f,V)=       I  I        Cmn(f)eJ2,fxmdXTeJ2,fYndvreJ2IIfYyT  ^ 

m  =  -MT    n  =  -NT 


where  cmn(f)  is  the  frequency  dependent  complex  weight  at  the  transmit  element  (m,n), 
and  dXx  and  dyr  are  the  interelement  spacings  in  the  X  and  Y  directions,  respectively. 
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( cross  range  ) 


(MRxNR) 


Fig.  2.3   Transmit  and  Receive  Array  Geometry 

The  last  complex  exponential  is  a  phase  factor  that  accounts  for  the  array  not  being  centered 
at  the  origin.    Also, 


MT  =  (  MT  -  1  )  /  2 


( 2.22a) 


and 


NT'  =  (  NT  -  1  )  /  2. 


(2.22b) 


It  should  be  noted  from  Eq.  (2.21)  that  a  single  omnidirectional  point  source  and  a  linear 
transmit  array  are  special  cases  of  the  planar  transmit  array.  This  is  also  what  we  would 
expect  from  examining  this  problem  from  a  purely  physical  point  of  view.  If  we  assume 
that  the  transmitted  electrical  signal  is  applied  to  the  array  elements  before  the  complex 
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weighting,  then  we  can  say  that  an  identical  electrical  signal  is  applied  to  all  the  elements  in 
the  array.  Accordingly,  Eq.  (2.21)  can  be  combined  with  Eq.  (2.19)  to  describe  the 
acoustic  input  to  the  ocean  medium. 

Finally,  what  is  meant  by  the  far-field  directivity  function?   It  can  be  shown  that  a  field 
point  at  a  range  r  is  considered  to  be  in  the  far  field  if 


r>^  (2.23) 


where  R  is  the  maximum  radial  extent  of  the  aperture.  For  our  planar  transmit  array,  this 
distance  is  given  by 

2  l}a 

R  =  [(Mrdxr)   +(NT,dYT)  J      .  (2.24) 

The  physical  significance  of  being  in  the  far  field  is  that  the  range  is  large  compared  to  a 
wavelength  and  the  separation  between  the  transducer  elements.  Consequently,  any 
second  or  higher  order  terms  in  the  expression  for  the  directivity  function  can  be  ignored 
and  the  mathematics  of  the  problem  is  simplified.  [Ref.  3:p.  38]  [Ref.  2:p.  170] 

F.   THE  PHASED  ARRAY 

Now  let  us  take  a  closer  look  at  the  frequency-dependent  complex  weighting  function 
Cmn  (0  of  the  far- field  directivity  function  given  by  Eq.  (2.21).  Since  this  weighting 
function  is  complex,  it  can  be  expressed  as 


Cmn(0  =  amn(f)ejemn(0  (2.25) 
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which  consists  of  an  amplitude  weighting  component  amn  (f)  and  a  phase  weighting  term 
®  mn  (0  .    Now  consider  that  this  phase  weighting  is  a  linear  phase  variation  of  the  form 

€>mn  (0  =  -2^  (  f  x  mdXT  +  f  y  nd yt  )  (2.26) 

where  m,  n,  dXT,  and  dyr  are  defined  as  for  Eq.  (2.21).  If  we  define  the  spatial 
frequency  terms  in  Eq.  (2.26)  as 

fx=u'A  (2.27a) 

and 

fy=v'A,  (2.27b) 

then,  substitute  Eqs.  (2.25)  and  (2.26)  into  Eq.  (2.21),  it  can  be  shown  that  the  resulting 
far-field  directivity  function  has  been  steered  in  the  direction  u  =  u'  and  v  =  v'.[Ref.  3:pp. 
132-134]. 

The  term  phased  array  means  nothing  more  than  an  array  which  employs  beamsteering 
using  phase  weighting.  Inspection  of  our  assumed  form  of  the  phase  weighting  term 
shows  that  it  is  both  a  separable  function  and  an  odd  function  of  the  indices  m  and  n. 
These  are  necessary  conditions  in  order  to  carry  out  the  beamsteering.  It  should  be  noted 
from  the  preceding  development  that  no  assumption  was  made  as  to  whether  the  amplitude 
term  was  either  separable  or  an  odd  function  of  the  indices  m  and  n.  For  the  purposes  of 
this  thesis,  we  will  only  be  considering  rectangular  amplitude  weighting,  which  is  a 
separable  function  and  leads  to  a  simplification  in  the  analysis. 


Since  the  phase  weighting  is  already  a  separable  function,  the  assumption  that  the 
amplitude  weighting  is  also  separable  results  in  the  entire  complex  weighting  term  given  by 
Eq.  (2.25)  being  separable.  Substitution  of  this  form  of  weighting  into  Eq.  (2.21)  results 
in  the  expression  for  the  far-field  directivity  function  being  separable  which  gives 

DT(f,v)  =  DT(f,fx)  DT(f,fY)  (2.28) 


where 


MT 

I 

m=-MT 


DT(f,fx)=     I     Aej2,tfXmd^e-j2'tfxmdCT,  (2.29) 


A  is  the  magnitude  of  the  amplitude  weighting,  and  a  similar  expression  exists  for 
DT(f,fy).  To  give  a  more  intuitive  feel  for  what  this  far-field  directivity  function  looks 
like,  we  can  express  Eq.  (2.29)  in  the  alternative  form 


^   ,cc    s      a  sin^fX  -fxJMTdxr] 

DT  (f  ,f  x )  =  A ; (2.30) 

sin[7t(fx  -fx)dxr] 


which  looks  similar  to  a  sine  function  when  viewed  in  direction  cosine  space.  The 
expression  for  DT  (f,f Y )  is  of  the  same  form  as  Eq.  (2.30)  except  that  it  is  multiplied  by  the 
last  complex  exponential  phase  factor  in  Eq.  (2.21)  that  accounts  for  the  array  not  being 
centered  at  the  origin. 


G.  RECEIVED  ACOUSTIC  FIELD 

We  next  need  to  calculate  the  acoustic  signal  at  some  spatial  location  r  =  (x,y,z)  after 
transmission  through  the  ocean  medium.  Applying  the  space-variant,  time-invariant 
impulse  response  of  the  ocean  medium  to  Eq.  (2.5)  yields 


,-=/:  / 


vm(M*)=  I  xM(t-T,r-r0)hM(x,r0;r)dxdr0   .  (2.31) 


If  we  express  the  input  signal  in  the  previous  expression  in  terms  of  its  frequency  and 
angular  spectrum  and  collect  terms  which  combine  with  the  impulse  response  to  produce 
the  transfer  function  of  the  ocean  medium,  we  obtain 

YM(f,r)  =  J       XM(f,v)HM(f,v;r)e"j27rv,rdv  (2.32) 


which  is  the  frequency  spectrum  of  the  output  acoustic  field  and  is  one  of  our  desired 
results. 

Next  we  need  to  consider  some  receive  aperture  to  convert  the  acoustic  energy  into 
electrical  energy.  For  the  purposes  of  this  thesis,  we  will  assume  that  the  receive  aperture, 
like  the  transmitter,  is  a  planar  array  that  is  characterized  by  a  complex  aperture  function 
AR(f,r).  Analogous  to  the  development  of  Eq.(2.16),  it  can  be  shown  that  the 
corresponding  output  electrical  signal  is  given  by 


y(t,r)  =  J        YM  (f ,r)A  R(f,r)ej2nft  df  .  (2.33) 


We  will  assume  that  our  receive  array  is  located  in  a  plane  parallel  to  the  XY  plane  and 
consists  of  MR  x  NR  (  odd  )  point  receivers  centered  at  (xR  ,yR  ,yR ).  We  will  further 
assume  that  the  array  elements  are  equally  spaced  in  the  X  and  Y  directions.  Due  to  the 
assumption  of  point  receivers,  the  aperture  function  can  be  expressed  as 


MR'       NR' 
AR(f,r)=     X        X    wmn(05[x-(xR4mdXR)]8[y-(yR+ndYR)]8[z-zR]     (2.34) 
m=-MR'  n=-NR' 


where  wmn  (f)  is  the  frequency  dependent  complex  weight  at  element  (m,n),  dxR  and  dyR 
are  the  interelement  spacings  in  the  X  and  Y  directions,  respectively, 

MR'  =  (  MR  -  1  )  /  2,  (2.35a) 


and 


NR'  =  (  NR  -  1  )  /  2.  (2.35b) 


For  this  thesis,  we  are  not  concerned  with  processing  the  received  electrical  signals. 
What  we  want  to  accomplish  is  the  ability  to  compute  the  complex  acoustic  field  at  the 
receive  array  as  a  function  of  frequency  and  spatial  coordinates,  Eq.  (2.32),  and  the  time- 
domain  electrical  signal  at  each  element  in  the  receive  array.  The  signal  processing  that 
will  be  done  on  the  computer  simulated  received  electrical  signals  will  be  via  auxiliary 
programs,  e.g.,  FFT  and  adaptive  beamformers.  Accordingly,  we  are  primarily  concerned 
with  formatting  our  outputs  in  a  form  acceptable  to  those  programs. 
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As  a  consequence  of  these  requirements,  and  in  an  attempt  to  minimize  computation 
time,  we  will  assume  the  simplest  form  of  the  aperture  function  where  the  complex 
weighting  function  wmn  (f)  is  unity.  Applying  this  simplifying  assumption  to  Eq.  (2.34) 
and  substituting  the  resulting  form  of  the  aperture  function  into  Eq.  (2.33)  yields 


;t,m,n,zR)=   I 


y(t,m,n,zR )  =  YM  (f,m,n,zR  )ej2*ftdf  (2.36) 


which  is  the  time  domain  electrical  signal  at  element  (m,n)  and  is  our  other  desired  result 

H.   OVERALL  SYSTEM  COMPLEX  FREQUENCY  RESPONSE 

If  we  combine  Eqs.  (2.19)  and  (2.32)  we  obtain 

YM(f,r)  =  X(f)H(f,r)  (2.37) 


where 


)  =   (      DT(f,v): 


H(f,r)  =  DT(f,v)HM(f,v;r)e"j27rV'rdv  (2.38) 


is  known  as  the  overall  system  complex  frequency  response  [Ref.  4:p.  5]  [Ref.  6:p.  1672]. 
Note  that  in  the  development  carried  out  in  this  chapter  we  have  not  made  any  assumptions 
about  the  ocean  medium  transfer  function.  The  only  assumption  used  in  the  derivation  of 
Eq.  (2.35)  was  that  an  identical  electrical  signal  was  applied  at  all  spatial  locations  of  the 
transmit  array  before  the  application  of  the  complex  weights. 
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Taking  the  inverse  time-domain  Fourier  transform  of  Eq.  (2.37),  and  using  the 
expansion  of  Eq.  (2.38)  it  can  be  shown  that 

YM  (V)  =(        f      X(f)DT(f,v)HM  (f,v;r)e-j27CV#r  ej27tfl  dv  df  (2.39) 


which  is  the  coupling  equation  that  relates  the  output  acoustic  signal  from  the  medium, 
y  M  (l'r) ,  to  the  frequency  spectrum  of  the  applied  elctrical  signal  X(f). 
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m.   THE  AXISYMMETRIC  CASE 

A.   OVERALL  SYSTEM  COMPLEX  FREQUENCY  RESPONSE  REVISITED 

Let  us  examine  the  expression  for  the  overall  system  complex  frequency  response, 
given  by  Eq.  (2.38),  that  was  developed  in  the  previous  chapter.  It  would  appear  from  an 
initial  inspection  that  Eq.  (2.38)  is  a  three-dimensional  integration  with  respect  to 
dv  =  df  x  df  y  df z  .  However,  it  should  be  noted  that  since  the  integrand  is  a  function  of 
the  frequency,  f,  as  well  as  the  spatial  frequencies,  v  =  ( fx.fy^z)'  tne  integral  can  be 
reduced  from  a  three-dimensional  to  a  two-dimensional  integral.  This  is  due  to  the  manner 
in  which  the  spatial  frequencies  given  by  Eq.  (2.7)  are  related  to  the  frequency  f  via 
direction  cosines.  The  property  of  direction  cosines  that  the  sum  of  their  squares  always 
equals  unity  means  that,  since  we  know  f,  we  can  express  one  of  the  direction  cosines  in 
terms  of  the  other  two.  Thus,  we  do  not  need  to  integrate  with  respect  to  all  three  spatial 
frequencies  in  order  to  span  direction  cosine  space.  Using  the  spatial  frequency  definitions 
of  Eq.  (2.7),  it  can  be  shown  that 

H(f,x,y ,z)  =    (        I        DT(f,f x  ,f Y  ,fz  )HM  (f,f x  ,fy,fZ  W.z) 

-j27t(f v  x  +  fYy  +  h z)  ,r     ,r 
x  eJ    vX        yj      z    dfxdfz       (3.1) 


where 


fY=±[(f/c0)2-(fx+f|)]      ,  (fx+f!)  ^(f/c0)2,  (3.2a) 

and 
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fy=  ^[(fl+f^-Cf/Co)2!     ,  (f|+fl)>(f/c0)2,  (3.2b) 


and  c0  is  the  speed  of  sound  at  a  particular  transmit  array  point  source  element.  The  plus 
(minus)  sign  in  Eq.  (3.2a)  is  chosen  whenever  y  >  y0  (y  <  yo)  and  corresponds  to  the 
field  point  being  deeper  (shallower)  than  the  source  point,  i.e.,  downward  (upward) 
traveling  waves.  The  minus  (plus)  sign  in  Eq.  (3.2b)  corresponds  to  the  plus  (minus) 
sign  of  Eq.  (3.2a)  in  order  to  generate  evanescent  waves. 

Note  that  we  have  not  made  any  assumptions  about  the  ocean  medium  transfer  function 
in  our  development  up  until  this  point.  However,  for  the  purposes  of  this  thesis,  we  are 
going  to  assume  that  the  acoustic  energy  of  the  transmitted  signal  is  contained  within  a 
sound  channel,  e.g.,  the  SOFAR  or  deep  sound  channel.  In  this  situation  the  sound 
propagates  between  upper  and  lower  turning  points  with  no  surface  and  bottom  boundary 
interaction.  We  therefore  expect  the  acoustic  field  trapped  within  this  channel  to  exhibit 
some  form  of  cylindrical  spreading.  If  this  assumption  is  correct,  then  cylindrical 
coordinates  should  provide  an  easy  reference  frame  against  which  to  describe  this  behavior. 

B.   THE  CYLINDRICAL  COORDINATE  SYSTEM 

The  cylindrical  coordinate  system  used  in  this  thesis  is  shown  in  Fig.  3.1,  where  the 
transformation  equations  that  relate  the  rectangular  and  cylindrical  coordinates  are  given  by 

x  =  rcoscj),  y  =  y,  andz  =  rsin<])  (3.3) 

I    2       2~~ 
where  r  =  v  x   +  z      is  the  polar  radial  distance. 

In  addition,  the  spatial  frequencies  must  be  transformed  using  the  following  set  of 

equations  (  see  Fig.  3.2  ): 
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X 

( cross  range  ) 


Fig.  3.1  The  Cylindrical  Coordinate  System. 


2itfz 


Fig  3.2  3-D  Wave  Propagation. 
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fx  =  frcos£,  fY  =  fy»  and  fz  =  fr  sin£  (3.4) 

V~2       2~ 
f  x  +  f  z    is  the  spatial  frequency  in  the  radial  direction.   As  a  result, 

dfxdfz->frdfrdC  (3.5) 

The  physical  significance  of  the  two  different  angles  in  Eqs.  (3.3)  and  (3.5)  is  that  <()  relates 
to  the  position  of  some  field  point  whereas  C  relates  to  the  direction  of  propagation,  and 
they  are  not  necessarily  the  same. 

Applying  these  definitions  to  Eq.  (3.1),  it  can  be  shown  that  the  overall  system 
complex  frequency  response  can  be  expressed  in  cylindrical  coordinates  as 


x  e-J27T(frr  cosC  sin*  +  fYy  +  frr  sinC  sine}))  f^^  (^6) 


H(f,r)=  |  I      DT(f,fr,C,fY)HM(f,fr,C,fY;r) 

'0        JQ> 


where 


and 


2      2    m        2  2 

fY=±[(f/c0)   -ft]      ,    fr2<(f/c0),  (3.7a) 


2  7l/2        2  2 

fy=  +j[fr  -(f/c0)  J      ,    fr>(f/c0),  (3.7b) 


and  Co  is  the  speed  of  sound  as  previously  defined. 
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Now  consider  the  specific  transmit  array  geometry  discussed  in  Chapter  2.  If  we 
substitute  the  far-field  directivity  function  for  this  array,  as  given  by  Eq.  (2.21),  into  the 
previous  equation,  and  perform  the  necessary  transformation  from  rectangular  to 
cylindrical  coordinates,  we  find  that 

(\2*      /•«         MT  NT 

H(i»  =  j  I  I  I      cmn(f)HM(f,fr,C,fY;r) 

•'O        •'O        m=-MT    n  =  -NT 
x  e"J27rfr(r  «*C  sin<|)  +  r  sin^  sin<}»  +  mdXTcosQ  e-j27tfY(y  -  yT  -  ndyrOf  df  ^         ~  g 

C.   THE  AXISYMMETRIC  OCEAN  MEDIUM  TRANSFER  FUNCTION 

In  this  thesis,  we  shall  assume  that  the  speed  of  sound  in  the  ocean  medium  is  a 
function  of  depth  only.  It  will  be  shown  in  the  following  chapter  that  under  this  condition 
the  ocean  medium  transfer  function  is  axisymmetric,  i.e.,  independent  of  the  azimuthal 
angle  (,  and,  as  a  result,  it  can  be  expressed  in  the  form 

HM(f,fr,C,fY;r,(t),y)  =  HM(f,fr,y0;y)  (3.9) 

where  y0  =yT  -  ndyj  is  tne  depth  of  a  point  source  in  the  nth  horizontal  row  of  the 
transmit  array.  Note  the  presence  of  the  terms  y  and  yG  in  Eq.  (3.9).  These  terms 
indicate  that  the  ocean  medium  transfer  function  is  dependent  on  the  speed  of  sound  at  the 
position  of  the  transmitter  as  well  as  the  receiver.  Since  the  assumed  form  of  our 
transmitter  is  a  planar  array,  we  need  to  consider  the  depth  and,  hence,  speed  of  sound  at 
each  element  in  the  array  when  calculating  the  acoustic  field  at  some  spatial  location  r.  This 
makes  physical  sense  when  considering  Eq.  (3.2),  which  shows  that  the  speed  of  sound  at 
the  transmit  element  affects  the  direction  of  propagation  and  the  threshold  between 
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propagating  and  evanescent  waves.     If  the  form  of  the  ocean  medium  transfer  function 
given  by  Eq.  (3.9)  is  substituted  into  Eq.  (3.8),  then  the  integrals  can  be  regrouped  to  give 

^oo        MT'  NT' 

H(f,D=/l  X  I      cran(f)HM(f,fr,y0;y)e-j2,tf^-5'°)fr 

/        e-j27cfr(r  cosC  sin<J>  +  r  sinC  sin<j>  +  mdXTcosO  ^^  „  .~ 

J n 


Next,  consider  the  integration  with  respect  to  £  in  the  preceding  equation.  Due  to 
our  placement  of  the  transmit  array  in  the  coordinate  system,  the  angle  <}>  is  the  horizontal 
polar  angle  (relative  to  the  positive  X  axis)  between  the  center  of  the  transmit  array  and  the 
field  point  of  interest.  Now  define  a  new  angle  $  m  that  is  the  horizontal  polar  angle 
between  a  specific  array  element  (m,n)  and  the  field  point,  as  shown  in  Fig.  3.3.  Note  that 
the  field  point  in  Fig.  3.3  is  in  three-dimensional  space,  but  that  it  is  only  the  (x,z) 
component  of  that  position  that  determines  <)).  Also  note  that  since  the  angle  is  in  the  XZ 
plane,  it  is  independent  of  the  vertical  position  y  and,  hence,  is  independent  of  the  array 
index  n.  Expressing  the  horizontal  polar  angle  <j)  in  this  integral  in  terms  of  the  new  angle 
<t>  m  and  collecting  terms,  results  in 

f       e-J27cfr(r  C0SC  sin*  +  r  sinC  sin<{>  +  mdXTcosO,c.  _    f         -j2jtfrrmcos(C  -  $m)tf       (3  1 1) 
•'o  •'O 


where 


2       2    m 
rm  =[(x-mdxr)   +z  ]  (3.12) 
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Fig  3.3  Graphical  Interpretation  of  <j)  m  (  Plan  View  ). 

is  the  polar  radial  distance  between  element  (m,n)  in  the  transmit  array  and  the  field  point, 
and  mdxr  is  the  horizontal  distance  between  element  (m,n)  and  the  center  of  the  array. 

Brekhovskikh  [Ref.  8:p.  76]  and  Officer  [Ref.  9:p.  126]  show  that  the  right-hand  side 
of  Eq.  (3.11)  is  in  the  form  of  the  integral  representation  of  the  zero-order  Bessel  function 
of  the  first  kind  where 


.( 


2ic 


-j2^frrrncos(C-0m)dC  =  27lJ()(Mrrm) 


(3.13) 


Substituting  Eqs.  (3.1 1)  and  (3.13)  into  Eq.  (3.10)  yields 


i-~      MT  NT 

H(f,r)  =  27t    I  X  X       cmn(0HM(f,fr,yo;y)J0(27cfrrm) 

•'O      m=-MT    n  =  -NT 

xe-j2rrfY(y-yn)frdfr 


(3.14) 
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which  is  the  fundamental  result  of  this  chapter.  Comparing  the  previous  equation  with  Eq. 
(3.1),  it  can  be  seen  that  our  expression  for  the  overall  system  complex  frequency  response 
has  been  reduced  from  a  double  to  a  single  integral.  This  was  based  on  the  assumption 
that  the  speed  of  sound  was  a  function  of  depth  only  and,  as  a  result,  the  corresponding 
ocean  medium  transfer  function  was  axisymmetric.  Since  numerical  integration  is  a  time 
consuming  process,  it  is  surmised  that  the  reduction  in  the  number  of  integrations  required 
in  the  computation  of  our  computer-generated  acoustic  output  should  result  in  the  saving  of 
a  significant  amount  of  computer  time. 

It  should  be  emphasized  that  Eq.  (3.14)  is  valid  for  any  axisymmetric  ocean  medium 
transfer  function.  Since  this  condition  is  satisfied  for  the  case  when  the  sound  speed  is  a 
function  of  depth  only,  Eq.  (3.14)  represents  the  general  form  of  the  overall  system 
complex  frequency  response  that  we  will  use  in  the  development  of  the  next  chapter. 
Remember  that  we  are  trying  to  develop  a  general,  modular,  pulse-propagation  model 
based  on  sound-speed  profiles  that  are  a  function  of  depth.  Thus  Eq.  (3.14),  when 
combined  with  the  coupling  equations  of  Chapter  2,  is  the  general  solution  for  this 
assumed  type  of  sound-speed  profile.  In  the  next  chapter  we  will  develop  a  specific 
solution  based  on  the  WKB  approximation. 

D.    AXIAL  SYMMETRY  OF  THE  TRANSMIT  ARRAY 

The  only  assumption  made  regarding  axial  symmetry  in  this  development  was  that  the 
ocean  medium  was  axisymmetric.  This  makes  physical  sense  when  considering  the  ocean 
medium  itself,  in  the  absence  of  any  source  of  acoustic  signal.  Since  the  speed  of  sound 
was  assumed  to  be  a  function  of  depth  only,  there  is  nothing  in  the  medium  that  should 
cause  any  azimuthal  <))  dependence  on  the  propagation  of  an  acoustic  signal.  Accordingly, 
the  ocean  medium  transfer  function  should  be  independent  of  the  azimuthal  angle  0. 
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The  assumption  of  axial  symmetry  may  at  first  glance  appear  to  be  invalid  when  we 
include  the  transmit  array  in  the  problem,  as  the  beam  pattern  or  far-field  directivity 
function  of  the  planar  array  is  not  generally  symmetrical  about  the  Y  axis.  However, 
consider  the  problem  from  the  following  perspective.  From  linear  systems  theory,  the  far- 
field  directivity  function  of  the  array  is  the  linear  superposition  of  the  far-field  directivity 
functions  of  the  individual  array  elements.  The  far-field  directivity  function  of  a  point 
source  is  unity,  i.e.,  it  is  omnidirectional.  The  far-field  directivity  function  of  a  complex 
weighted  point  source  is  also  omnidirectional,  since  the  complex  weight  factor  is  merely  a 
scaling  term.  Equations  (2.20),  (2.21)  and  (2.29)  demonstrate  that  the  far-field  directivity 
function  of  the  transmit  array  is  a  result  of  the  way  these  omnidirectional  beam  patterns  are 
combined  rather  than  by  altering  the  beam  pattern  of  the  individual  point  sources.  The 
effect  of  the  complex  weighting  is  to  control  the  way  these  individual  far-field  directivity 
functions  are  combined. 

Examination  of  Eq.  (2.21)  also  shows  that  there  is  an  additional  phase  term  that  is 
associated  with  each  array  element  that  accounts  for  that  element's  displacement  from  the 
origin  of  the  coordinate  axis  system.  Thus,  we  can  assume  that  the  far-field  directivity 
function  of  the  transmit  array  is  the  linear  combination  of  complex  weighted  axisymmetric 
point  sources.  Accordingly,  the  overall  system  complex  frequency  response  solution 
given  by  Eq.  (3.14)  can  be  thought  of  as  applying  to  each  individual  array  element,  with 
the  solution  for  the  array  being  a  linear  combination  of  these  results.  This  is  nothing  more 
than  rewriting  Eq.  (3.14)  with  the  summations  outside  of  the  integral.  Since  the 
integration  and  summations  are  with  respect  to  different  quantities,  regardless  of  whether 
the  integration  is  inside  the  summations  or  outside,  the  value  of  the  overall  system  complex 
frequency  response  will  be  the  same. 


IV.   THE  OCEAN  MEDIUM  TRANSFER  FUNCTION 

A.   THE  HELMHOLTZ  WAVE  EQUATION 

Examining  the  various  coupling  equations  developed  in  Chapter  2,  we  see  that  they 
already  incorporate  an  input  acoustic  signal  xM  (t,r),  with  corresponding  frequency  and 
angular  spectrum  XM  (f,v).  This  means  that  the  transfer  function  that  characterizes  the 
ocean  medium,  HM(f,v;r),  is  independent  of  the  input  as  we  would  expect  from  linear 
systems  theory.  Accordingly,  when  trying  to  find  a  transfer  function  to  describe  the 
ocean  medium  we  only  need  to  solve  the  wave  equation  given  by 


V2(p(t,r)-^— ~cp(t,r)-0  (4.1) 

c  (r)  at2 


instead  of  the  linear  acoustic  wave  equation  given  by  Eq.  (2.1). 

The  first  assumption  that  we  will  make  in  finding  a  solution  to  Eq.  (4.1)  is  that  the 
velocity  potential,  cp^1"),  has  a  time-harmonic  dependence,  that  is, 


cp(t,r)  =  (p(r)ej27rft.  (4.2) 


If  this  form  of  the  velocity  potential  is  then  substituted  into  Eq.  (4.1)  and  the  time-harmonic 
dependence  is  factored  out,  the  result  is  the  time-independent  Helmholtz  wave  equation 
given  by 


V2(p(r)  +  k2(r)(p(r)  =  0  (4.3) 
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where  k(r)  =  2rcf  /c(r)  is  the  wave  number  and  c(r)  is  the  speed  of  sound  in  the  fluid 
medium  at  spatial  location  r. 

As  discussed  in  the  previous  chapter,  we  will  be  assuming  for  the  purposes  of  this 
thesis  that  the  speed  of  sound  is  a  function  of  only  one  spatial  variable,  the  depth  y,  i.e., 
c(r)  =  c(y).  While  this  assumption  is  not  in  general  true,  it  is  nevertheless  adequate 
enough  to  describe  the  ocean  medium  for  most  circumstances  where  we  are  not  considering 
propagation  ranges  longer  than  say,  20  km,  or  propagation  across  significant 
oceanographic  perturbations  such  as  fronts  or  eddies.  We  can  also  consider  that  the 
sound- speed  profile  in  the  region  of  the  deep  sound  channel  is  less  variable  than  for 
shallow  sound  channels  or  surface  ducts.  Thus,  the  assumption  that  the  sound-speed 
profile  is  not  dependent  on  geographical  position  is  adequate  for  our  model,  within  the 
limitations  discussed. 

This  assumption  allows  us  to  make  the  substitution  k(r)  =  k(y)  in  Eq.  (4.3)  which 
greatly  simplifies  the  analysis.  Since  k(r)  is  now  a  function  of  only  one  variable,  we  can 
use  the  method  of  separation  of  variables  to  solve  for  the  acoustic  velocity  potential. 
Applying  this  solution  technique  in  cylindrical  coordinates,  we  assume  a  solution  for  (p(r) 
of  the  form 

<p(r)  =  R(r)0>«j>)Y(y)  (4.4) 

where  r  =  (r,0,y)  is  the  spatial  location  in  cylindrical  coordinates  (see  Fig.  3.1)  and  the 
Laplacian  appearing  in  Eqs.  (4.1)  and  (4.3)  is  defined  as  being 


V2=il  +  lA+>  A_  +  ^.  (4.5) 

dr       r  dr      r    dty       dy 
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B.  AZIMUTHAL  ANGLE  AND  RANGE  TERM  SOLUTION 

Since  the  speed  of  sound  is  a  function  of  depth  only,  there  is  nothing  in  the  medium 
that  should  cause  any  azimuthal,  i.e.,  <J>,  dependence  on  the  propagation  of  an  acoustic 
signal.  Accordingly,  the  solution  for  the  acoustic  velocity  potential  should  be  independent 
of  the  azimuthal  angle  <|>.  This  is  referred  to  as  the  ocean  medium  having  axial  symmetry. 
Applying  this  assumption,  Eq.  (4.4)  reduces  to 

q>(r)  =  R(r)Y(y)  (4.6) 

and  the  partial  differential  with  respect  to  azimuthal  angle  in  the  Laplacian  of  Eq.  (4.5) 
drops  out. 

If  we  substitute  Eq.  (4.4)  into  Eq.  (4.3)  and  assume  that  k(r)  =  k(y),  the  method  of 
separation  of  variables  yields 


9_R«+IiRfi)+krR(r)  =  0  (4.7) 

9r2         r     dr 


where  the  wave  number  in  the  radial  direction,  kr ,  has  been  used  as  a  separation  constant. 
Making  the  substitution  R(r)  =  g(  krr  )  into  the  previous  equation,  it  can  be  shown  that  Eq. 
(4.7)  reduces  to  Bessel's  differential  equation  of  order  zero.  The  solution  to  this  equation 
can  be  expressed  in  either  of  the  forms 

R(r)  =  A  J0  (krr)  +  B  Y0  (krr)  (4.8a) 
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R(r)  =  A  H^(krr)  +  B  Hf,2)(krr)  (4.8b) 


where  Jo  and  Y0  are  the  zero  order  Bessel  functions  of  the  first  and  second  kind, 
respectively,  with  Yq  also  being  known  as  the  zero  order  Neumann  function.  The 
functions  H0  and  H0  are  zero  order  Hankel  functions  of  the  first  and  second  kind, 
respectively,  or  Bessel  functions  of  the  third  kind.  Since  Eqs.  (4.8a)  and  (4.8b)  are 
alternative  expressions  describing  the  same  quantity,  the  Hankel  functions  should  be  able  to 
be  expressed  in  terms  of  the  Bessel  functions  of  the  first  and  second  kind  and  visa-versa. 
This  is,  in  fact,  the  case  and  we  will  make  use  of  one  of  these  relationships  at  a  later  stage 
of  the  analysis. 

Since  there  are  two  alternative  representations  of  the  range  term  R(r),  we  would  expect 
that  one  representation  is  more  useful  than  the  other  for  the  situation  that  we  wish  to  model. 
If  we  consider  the  approximations  for  large  argument  for  the  two  Hankel  functions 


,,(1),.     ,  I     2        j(kr  -n/4) 

H0Urr)^Y^7e  (49a) 


and 


-,(2),.      ,  /     2        -j(kr-ic/4) 

HoUrr)-^V^rTe  .  (4-9b) 


then  because  of  our  assumed  form  of  the  time-harmonic  dependence  of  the  velocity 
potential  given  by  Eq.  (4.2),  Eq.  (4.9a)  represents  an  incoming  wave  while  Eq.  (4.9b) 
represents  an  outgoing  wave.  Since  we  are  assuming  that  our  transmitted  acoustic  signal 
is  traveling  in  a  channel  with  no  boundary  interaction,  we  will  have  only  outgoing  waves 


35 


from  the  source  and  no  incoming  waves  due  to  reflection.  It  then  appears  that  Eq.  (4.8b) 
with  the  constant  A  set  to  zero,  is  an  appropriate  solution  for  our  range  term.  We  will 
therefore  assume  that  our  range  term  solution  is  of  the  form 


R(r)  =  BrH^2)(krr)  (4.10) 


where  Br  is  a  constant.  By  comparison,  if  we  look  at  the  approximations  for  large 
argument  of  the  Bessel  functions  of  the  first  and  second  kind,  we  find  that  they  are 
comprised  of  both  incoming  and  outgoing  wave  components  and  so  are  not  well  suited  to 
our  problem. 

Upon  substitution  of  Eq.  (4.10)  into  Eq.  (4.7)  however,  we  find  that  our  assumed 

form  of  the  solution  for  the  range  term  is  a  solution  everywhere  except  at  the  origin,  i.e., 

(2) 
r  =  0,  since  Hq   "blows  up"  at  this  point  and  equals   -j  ©o  .     Using  a  volume  integral 

approach,  Gauss's  theorem,  and  the  asymptotic  expression  for  the  Hankel  function  given 

by 


lim     H^2)(krr)=^ln(krr),  (4.11) 

r->0  nT 


it  can  be  shown  that  Eq.  (4.10)  is  a  solution  to 


[V2+kr2]R(r)  =  ^L5o-)  (4.12) 


for  all  r,  where  Eq.  (4.12)  is  the  same  as  Eq.  (4.7)  with  the  addition  of  a  forcing  function 
or  source  term.    This  makes  physical  sense  when  we  remember  that  Eq.  (4.10)  describes 
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outgoing  waves  in  the  radial  direction.  Since  these  outgoing  waves  have  to  come  from 
somewhere  and  there  are  no  incoming  waves,  there  must  be  a  source  of  some  description, 
that  generates  the  waves.  Thus  Eq.  (4.10)  is  a  solution  to  Eq.  (4.12)  that  satisfies  the 
boundary  condition  of  a  point  source  at  the  origin.  The  choice  of  a  point  source  as  the 
forcing  function  is  appropriate  since  we  have  assumed  that  our  transmit  aperture  is  a  planar 
array  of  point  sources. 

C.    THE  WKB  APPROXIMATION 

For  the  depth  dependent  term  Y(y)  of  our  solution  for  the  velocity  potential,  we  will 
obtain  an  approximate  solution  based  on  the  geometrical  optics  approximation,  otherwise 
known  as  the  WKB  (Wentzel,  Kramers,  and  Brillouin)  approximation.  If  we  substitute 
Eq.  (4.4)  into  Eq.  (4.3)  and  assume  that  k(r)  =  k(y),  the  method  of  separation  of  variables 
yields 


[^+k2Y(y)l  Y(y)  =  0  (4-13) 

dy 


where  the  wave  number  in  the  vertical  (or  depth)  direction,  kY  (y),  is  calculated  using  the 
relationship 


kY(y)  =  k2(y)-k2.  (4.14) 


If  the  vertical  wave  number  was  a  constant,  i.e.,  kY(y)  =  kY,  meaning  that  it  was  no 
longer  a  function  of  depth,  then  Eq.  (4. 13)  would  have  the  exact  solution 

Y(y)  =  Ae-jkYy+BejkYy.  (4.15) 


If  we  now  assume  that  the  local  medium  is  homogeneous,  that  is  the  variations  of  the 
properties  of  the  medium  per  wavelength  are  small,  then  a  solution  similar  in  form  to  Eq. 
(4.15)  would  be  a  good  approximation.  The  WKB  method  assumes  that  the  solution  to 
Eq.  (4.13)  is  of  the  form 

Y(y)  =  AY(y)ejeY(y)  (4.16) 

where  AY  (y)  and  0  Y  (y)  arc  rca^  amplitude  and  phase  functions  respectively.  Substituting 
this  assumed  form  of  our  solution  into  Eq.  (4.13),  and  making  the  simplifying  assumption 
that 


|AY"(y)/AY(y)|«[0Y'(y)]2  (4.17) 


yields 


v, ,         l rA     -jfykY(OdC         j/"ykY(Qdi;  -i 


which  is  the  approximate  WKB  solution  [Ref.3:pp.  208-220]  [Ref.  8:pp.  129-134].  It 
should  be  noted  that  the  only  assumption  used  that  leads  to  this  not  being  an  exact  solution 
is  that  given  by  Eq.  (4.17).  This  method  also  assumes  that  there  is  no  reflection  in  a 
horizontally  stratified  medium. 

The  LHS  of  Eq.  (4.17)  is  the  relative  differential  of  the  time  rate  of  change  of 
amplitude  and  is  a  measure  of  the  rate  of  change  in  the  properties  of  the  medium.  If  the 
medium  was  homogeneous,  then  this  term  would  be  zero  since  the  characteristics  of  the 
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medium  would  not  be  varying.  Therefore  this  term  is  a  measure  of  the  inhomogeneity  of 
the  medium.  The  RHS  of  Eq.  (4.17)  is  inversely  proportional  to  wavelength  since  for  a 
given  frequency  the  rate  of  change  of  phase  increases  with  decreasing  wavelength. 
Equation  (4.17)  is  therefore  a  reasonable  simplifying  assumption  since  it  is  nothing  more 
than  a  mathematical  formulation  of  the  assumption  made  earlier  that  the  variations  of  the 
properties  of  the  medium  are  small  compared  to  a  wavelength. 

The  integrals  in  Eq.  (4.18)  give  the  phase  change  as  the  waves  propagate  between 
depths  y  and  yG,  while  the  factor  in  front  is  to  preserve  the  conservation  of  energy.  Due  to 
the  assumed  form  of  the  time-harmonic  dependence  of  the  velocity  potential  given  in  Eq. 
(4.2),  the  first  and  second  terms  in  Eq.  (4.18)  correspond  to  waves  traveling  in  the  positive 
(down)  and  negative  (up)  Y  directions,  respectively.  It  should  be  noted  that  the  WKB 
solution  is  then  the  superposition  of  two  waves  traveling  without  interaction  in  opposite 
directions.  Also  note  that  the  approximate  solution  is  not  valid  as  ky  (y)  approaches  zero 
since  the  solution  approaches  infinity.  The  point  at  which  this  occurs  is  referred  to  as  a 
turning  point. 

D.   THE  VELOCITY  POTENTIAL  SOLUTION 

Combining  Eqs.  (4.6),  (4.10)  and  (4.18),  we  obtain  a  solution  for  the  velocity  potential 
of  the  form 


AU(o\krr)     -J   f   MOdC 
— ==—  e    •'v. 


9(r)-  #yly) 


(4.19) 


where 


1  /2 
kY(y)  =  ±27t[(f/c(y))2-f?]         ,    (fr)2<(f/c(y))2,  (4.20a) 
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and 


2  2     1/2  2  2 

ky  (y)  =  +j2tc[  fr  -  (  f /c(y)  )Z  ]         ,    (fr  )   >  (  f /c(y) )   ,  (4.20b) 


which  is  our  separable  function  solution  to  the  wave  equation.  The  constant  A  in  Eq. 
(4. 19)  combines  the  constants  Ay  and  B  r  of  the  preceding  equations.  Note  that  we  have 
replaced  the  approximation  symbol  with  the  equality  sign  in  Eq.  (4.19).  This  is  only  a 
matter  of  notation,  and  we  still  recognize  that  the  solution  is  only  an  approximation. 

Since  the  linear  wave  equation  is  in  terms  of  partial  differentials  with  respect  to  y  and  r, 
if  we  multiply  Eq.  (4.19)  by  a  function  of  kr,  then  the  result  is  also  a  solution  to  the  wave 
equation  since  kr  is  independent  of  y  and  r  [Ref.  9:p.  125].  Similarly,  any  summation  of 
such  solutions  will  also  be  a  solution  to  the  wave  equation.  Combining  these  concepts  into 
a  generalized  integral  form  results  in 

f°°   A(kr)H?\krr)     -J  f  MOdC 

where  A(kr)  is  the  term  that  incorporates  the  constant  A  of  Eq.  (4.19)  and  the  arbitrary 
function  of  kr  as  discussed.  Equation  (4.21)  is  therefore  our  solution  of  the  linear  wave 
equation  based  on  the  WKB  solution  for  an  inhomogeneous  medium.  If  Eq.  (4.21)  is 
correct,  then  for  the  case  of  a  homogeneous  medium  it  should  collapse  into  the  free-space 
Green's  function  for  an  isospeed  medium.  As  well  as  validating  the  form  of  the  equation, 
it  should  also  provide  us  with  an  expression  for  the  term  A(kr). 
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E.  THE  FREE-SPACE  GREEN'S  FUNCTION 

The  free-space  problem  for  an  isospeed  medium  is  one  in  which  there  are  no 
boundaries  to  consider  and  the  medium  is  considered  to  be  homogeneous.  Since  one  of 
the  assumptions  made  in  the  development  of  Eq.  (4.21)  was  that  there  was  no  surface  or 
bottom  boundary  interaction,  all  that  we  have  to  do  to  find  a  solution  for  the  free-space, 
isospeed  problem  is  to  apply  the  conditions  for  a  homogeneous  medium  to  Eq.  (4.21). 
For  the  homogeneous  medium,  kY  (y)  =  kY  (y0)  is  a  constant.  For  this  case,  there  is  a 
closed  form  solution  of  the  integral  appearing  in  the  exponent  of  Eq.  (4.21),  which  results 
in  a  solution  of  the  form 

m=    r~    A(k)H0\r)  e-jkY<y0)<y-yo)dki  (4  22) 

•'-co       ^|kY(y0)| 

It  should  be  noted  that  for  the  case  of  a  homogeneous  medium,  the  simplifying  assumption 
given  by  Eq.  (4.17)  is  unnecessary,  and  Eq.  (4.22)  is  an  exact  solution. 

We  will  now  specifically  solve  the  time-independent  Helmholtz  equation  of  the 
form  of  Eq.  (4.3)  for  the  free-space,  isospeed  problem  and  compare  the  result  with  Eq. 
(4.22).  In  the  free  space  problem  a  wave  will  radiate  in  an  outwards  direction, 
experiencing  spherical  spreading,  and  be  free  of  boundary  interactions.  The  homogeneous 
medium  condition  means  that  k(r)  =  k  is  a  constant.  Since  the  wave  will  spread 
spherically,  it  makes  sense  to  use  the  spherical  coordinate  system  to  describe  the  solution. 
It  can  be  shown  that 


A  "J'kr 
cp(r)  =  -^—  (4.23) 


is  a  solution  to  Eq.  (4.3),  under  the  conditions  described,  everywhere  except  at  r  =  0. 
Note  that  in  Eq.  (4.23)  r  is  a  scalar  quantity  rather  than  a  vector.  Equation  (4.23), 
however,  does  satisfy 


V  2<p(r)  +  k2<p(r)  =  -47iA8(r)  (4.24) 


for  all  r.  This  is  consistent  with  there  having  to  be  a  source  at  r  =  0  to  generate  the  wave. 
The  choice  of  a  point  source  as  the  source  is  a  convenient  one  since  in  Chapter  2  we 
assumed  that  our  transmit  aperture  was  an  array  of  point  sources.  If  we  think  of  this 
forcing  function  or  source  term  as  having  unit  amplitude,  we  set  the  constant  in  Eq.  (4.23) 
equal  to  (-1/4tu).  Applying  this  amplitude  term  and  converting  to  cylindrical  coordinates 
we  can  show  that 


is  a  solution  to 


for  all  r  where 


-jkR 
<P(r)  =  ^  (4.25) 


V2cp(r)  +  k2cp(r)  =  ^  5(y-y0 )  (4.26) 


2  2    1/2 

R  =  [r    +(y-y0)   ]      .  (4.27) 
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and  r  is  the  polar  radial  distance.  Equations  (4.23)  and  (4.25)  are  the  free-space  Green's 
functions  for  an  isospeed  medium  in  spherical  and  cylindrical  coordinates,  respectively. 
They  represent  the  solution  to  the  time-independent  Helmholtz  equation  for  a  homogeneous 
medium  with  no  boundary  interactions,  and  matching  the  boundary  condition  for  a  point 
source. 

Officer  [Ref.  9:p.  126]  shows  that  the  free-space  Green's  function  in  cylindrical 
coordinates  can  be  expressed  as  an  integral  of  the  zero-order  Bessel  function  of  the  first 
kind,  having  the  form 


-jkR        r  oo        . 
V  =  J      Rfe'oiMe**-^*,.  (4.28) 


If  we  now  combine  Eqs.  (4.28)  and  (4.25)  we  obtain 

*r)  =  f~  ^A_  Jo(M)  e^oXy-yo)^  (429) 

which  is  the  solution  found  by  solving  the  linear  wave  equation  specifically  for  the  free- 
space,  isospeed  problem. 

F.   DETERMINATION  OF  THE  OCEAN  MEDIUM  TRANSFER  FUNCTION 

Returning  to  our  solution  for  the  free-space  problem  based  on  the  WKB  approximation 
given  by  Eq.  (4.22),  we  see  that  it  is  in  a  similar  form  to  Eq.  (4.29)  but  is  in  terms  of  a 
Hankel  function  rather  than  a  Bessel  function.    Based  on  observations  made  earlier  in  this 

chapter,  however,  we  would  expect  to  be  able  to  express  Eq.  (4.22)  in  terms  of  a  Bessel 

2  1 

function.    This  is  in  fact  the  case,  and  using  the  relationships  Ho(-krr)  =  -Ho(krr)  and 
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J0(krr)  =  i[Hj(krr)  +  Ho(krr)]from  Brekhovskikh  [Ref.  8:pp.  76-77]  it  can  be 
shown  that 

=  p2^o^c,kY(yo)(y.yo) 

J0      #v(yo)| 


Our  two  solutions,  given  by  Eqs.  (4.29)  and  (4.30),  are  now  in  a  comparable  form  and 
we  can  determine  the  value  of  the  constant  A(kr )  to  be 


«^ 


Note  the  fact  that  the  two  terms  involving  the  wave  number  in  the  vertical  or  depth 
direction,  kY(y0),  in  Eq.  (4.31)  are  not  partially  cancelled.  This  is  due  to  the  fact  that 
their  value  may  be  either  purely  real  (propagating  wave)  or  purely  imaginary  (evanescent 
wave)  caused  by  the  infinite  range  of  kr.  As  a  consequence,  Eq.  (4.31)  is  the  most 
concise  way  of  expressing  this  quantity. 

The  fact  that  we  were  able  to  get  our  solution  from  the  WKB  approximation,  Eq. 
(4.30),  into  the  same  form  as  the  Green's  function,  Eq.  (4.29),  for  the  free-space,isospeed 
problem  suggests  that  our  form  of  the  WKB  solution  is  correct.  Combining  Eqs.  (4.31), 
(4.21),  and  changing  the  variable  of  integration  from  krto  fr  it  can  be  shown  that 


cp(r): 


f~      J7tVkY(y0)  -jf      kY(C)dC 

I       v    t     \   /IT— ^=fJ0(M)e     JVo  frdfr>  (4-32) 

Jn    kY(y0)VkY(y) 
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which  is  our  solution  based  on  the  WKB  approximation  for  a  sound-speed  profile  that  is  a 
function  of  depth  only.  This  is  one  of  the  fundamental  results  of  this  chapter.  Now  that 
this  quantity  is  known,  we  can  calculate  the  ocean  medium  transfer  function  based  on  the 
WKB  approximation. 

The  simplest  way  to  calculate  the  ocean  medium  transfer  function  is  to  compare  Eq. 
(4.32)  with  the  overall  system  complex  frequency  response  of  Eq.  (3.14)  calculated  for  a 
single  point  source  of  unit  amplitude.  Under  these  circumstances  MT,  NT,  and  cmn(f) 
equal  unity  and  Eq.  (3.14)  reduces  to 

H(f,r)  =  27t    f     HM(f,fr,yo;y)J0(27tfrrm)e-j27rfY(yo)(y-yn)frdfr.  (4.33) 

The  reason  that  we  are  interested  in  Eq.  (4.33)  is  that  these  are  the  very  same  conditions 
that  led  to  the  solution  given  by  Eq.  (4.32).  Consequently,  Eqs.  (4.33)  and  (4.32)  should 
be  the  same.  Any  differences  between  these  two  equations  must  therefore  be  accounted 
for  by  the  ocean  medium  transfer  function  as  this  is  the  only  unknown  quantity  in  the 
equations.    It  can  therefore  be  shown  that 


HM(f.fr.yo;y)=  ^p^4>Me^vW-»J         (4.34) 
2kY(yo)v|kY(y)| 


where  kY(y)  is  as  defined  by  Eqs.  (4.20a)  and  (4.20b).  Equation  (4.34),  then,  is  our 
ocean  medium  transfer  function  based  on  the  WKB  approximation  and  is  the  fundamental 
result  of  this  chapter. 
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V.  THE  SIGNAL  GENERATOR 

A.  THE  COMPLEX  ENVELOPE 

Consider  a  transmitted  electrical  signal  of  the  following  form  : 

x(t)  =  a(t)cos(  27tfct  +  9(t) )  (5.1) 

This  is  an  amplitude  and  angle  modulated  carrier  where  we  will  assume  that  a(t)  and  @(t) 
are  arbitrary,  real  baseband  amplitude  and  angle  modulating  signals,  respectively.  This 
implies  that  both  a(t)  and  0(t)  are  bandlimited  to  the  region  -W  <  f  <  W.  These  signal 
components,  a(t)  and  0(t),  are  commonly  known  as  the  natural  envelope  and  phase  of  the 
signal  x(t),  respectively.  As  a  consequence,  x(t)  is  a  real  bandpass  or  bandlimited 
signal.  This  means  that  the  amplitude  spectrum  of  the  signal  is  zero  outside  a  frequency 
band  of  width  2W  Hz  centered  about  +/-  fc  Hz,  where  fc  >  W  Hz  and  fc  is  referred  to  as 
the  center  or  carrier  frequency. 

Using  Eq.  (5.1)  for  the  model  of  x(t),  we  will  make  use  of  the  concept  of  representing 
a  real  bandpass  signal  in  terms  of  its  complex  envelope  [Ref.  3:pp.l76-188]  [Ref.  10:pp. 
74-90].    The  complex  envelope  of  x(t),  denoted  by  x(t),  is  defined  as 


x(t)  =  [x(t)+jx(t)]e"j27rfct  (5.2) 


where  x(t)  is  the  Hilbert  transform  of  x(t).  If  x(t)  is  real,  then  its  Fourier  transform 
exhibits  Hermitian  symmetry.  Consequently,  the  frequency  spectrum  along  the  negative 
frequency  axis  is  completely  determined  by  the  spectrum  for  positive  frequencies  (and  visa- 
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versa).     Thus,  if  we  remove  the  frequency  spectrum  along  the  negative  (or  positive) 
frequency  axis,  we  will  not  loose  any  of  the  information  contained  in  the  signal. 
The  Hilbert  transform  of  x(t)  is  defined  by 

x(t)=l  f    : 


■?f  — * 

■      I  t-T 


which  is  the  convolution  of  x(t)  with  the  time  function  1/7U.  This  has  the  effect  of 
producing  a  -/+  90°  phase  shift  for  all  positive/negative  frequencies  of  the  input  signal. 
The  Hilbert  transform  applies  to  any  signal  that  is  Fourier  transformable,  and  by  looking  at 
Eq.  (5.3)  in  the  frequency  domain,  we  obtain  the  analogous  definition 

x(t)  =  IFT  {  -jsgn(f)  X(f)  }  (5.4) 

where  sgn(f)  is  the  signum  function  and  IFT{-}  is  the  inverse  Fourier  transform.  The 
effect  of  adding  a  real  signal  to  its  Hilbert  transform  is  to  remove  the  negative  frequency 
components,  producing  a  one-sided  frequency  spectrum  which  is  known  as  the 
pre-envelope  x+  (t).  [Ref.  10:pp.  74-75]  [Ref.lhpp.  67-72] 

This  one-sided  spectrum  is  then  translated  down  in  frequency  so  as  to  be  centered  at  0 
Hz  to  obtain  the  complex  envelope.  Thus,  the  complex  envelope  of  a  real  bandpass  signal 
is  typically  a  complex  signal  with  a  lowpass  or  baseband  frequency  spectrum  that  has  the 
same  information  content  as  the  original  bandpass  signal. 

Applying  the  definition  of  the  Hilbert  transform  to  the  specific  case  of  the  real  bandpass 
signal  x(t)  in  the  form  of  Eq.  (5.1)  results  in 
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x(t)  =  xc(t)sin(27tfct)  +  xs(t)cos(27rfct)  (5.5a) 

where 

xc(t)  =  a(Ocos0(t),  (5.5b) 

xs(t)  =  a(t)sin0(t),  (5.5c) 

and  xc  (t)  and  xs  (t)  are  known  as  the  in-phase  and  quadrature  components  of  the  bandpass 
signal,  respectively.  It  should  be  reiterated  here  that  one  of  the  underlying  assumptions  in 
this  treatment  of  Eq.  (5.1)  is  that  x(t)  is  a  bandpass  signal.  This  in  turn  requires  xc  (t)  and 
xs  (t)  to  be  bandlimited  to  the  interval  -W  <  f  <  W  where  fc  >  W.  Equation  (5.5a)  is 
sometimes  referred  to  as  the  canonical  form  of  x(t).  [Ref.  3:p.l83]  [Ref.  10:p.  85] 
By  combining  Eqs.  (5.1),  (5.2)  and  (5.5),  it  can  be  shown  that 


x(t)  =  a(t)ej6(t)  (5.6) 


is  the  complex  envelope  of  the  assumed  form  of  the  transmitted  electrical  signal  given  by 
Eq.  (5.1).  Given  our  initial  assumptions,  this  means  that  the  amplitude  spectrum  of  the 
complex  envelope  described  by  Eq.  (5.6)  will  be  baseband  centered  around  0  Hz  with  a 
bandwidth  of  W  Hz. 

The  form  of  Eq.  (5.6)  demonstrates  the  advantages  of  complex  envelope  notation. 
First,  it  allows  for  the  simple  representation  of  amplitude  and  angle  modulation.  Second, 
for  our  assumed  form  of  x(t)  given  by  Eq.  (5.1),  which  was  a  modulated  carrier,  the 
complex  envelope  is  independent   of  the  carrier  frequency.     The  significance  of  this 
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independence  from  the  carrier  becomes  apparent  when  dealing  with  the  discrete 
representation  of  x(t),  required  for  signal  processing  by  digital  computers.  Since  the 
information  content  of  the  signal  x(t)  is  completely  represented  by  its  complex  envelope 
x(t),  we  can  sample  the  baseband  complex  envelope  rather  than  the  original  bandpass 
signal.  As  the  maximum  frequency  component  present  in  x(t)  is  W  Hz  rather  than 
(fc  +  W)  Hz  as  in  x(t),  a  lower  sampling  rate  is  required  to  satisfy  Nyquist's  criterion. 
Consequently,  less  time  samples  or  data  points  will  be  required  to  represent  the  signal 
resulting  in  a  reduction  of  computation  time. 

It  can  be  shown  that  both  the  original  bandpass  signal  x(t)  and  its  Fourier  transform 
X(f)  can  be  expressed  in  terms  of  the  complex  envelope  using 


x(t)=Re[x(t)ej27rfcl]  (5.7) 


and 


X(f)=y[X(f-fc)+X(-(f-fc))].  (5.8) 


From  the  previous  equations  it  can  be  seen  that  it  is  a  straight  forward  matter  to  move 
between  the  complex  envelope  and  common  signal  representations  of  x(t)  in  both  the  time 
and  frequency  domains.  Thus  the  idea  of  working  with  the  complex  envelope  is  not  at  the 
expense  of  any  substantial  computational  penalty  in  going  from  one  representation  of  the 
signal  to  the  other.  [Ref.  lO.pp.  85-86] 
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B.  THE  CW  AND  LFM  PULSE 

For  the  purposes  of  this  thesis,  two  classic  waveforms  used  in  underwater  acoustics 
will  be  modelled.  They  are  the  CW  (Continuous  Wave)  and  LFM  (Linear  Frequency 
Modulated)  pulses,  both  of  which  are  real  bandpass  signals  that  can  easily  be  expressed  in 
the  general  form  of  Eq.  (5.1). 

A  CW  waveform  is  simply  one  in  which  there  is  no  angle  modulation.  Thus,  the 
expressions  for  the  bandpass  CW  signal  and  its  complex  envelope  can  easily  be  found  by 
setting  0(t)  equal  to  zero  in  Eqs.  (5.1)  and  (5.6)  yielding 

Xcw(t)=a(t)cos(27tfct),  (5.9) 


and 


Xcw(t)  =  a(t).  (5.10) 

Frequency  Modulation  is  one  of  the  most  commonly  used  types  of  angle  modulation 
and  is  where  the  instantaneous  frequency,  fi  (t),  is  varied  linearly  with  some  arbitrary 
baseband  signal.  The  instantaneous  frequency  is  related  to  the  time  derivative  of  the 
instantaneous  phase  and,  for  the  real  bandpass  signal  of  the  form  of  Eq.  (5.1),  it  can  be 
shown  that 


J_  d(8(t)) 
dt 


where  fc  is  the  carrier  frequency.  [Ref.  10:pp.  180-183] 
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(5.11) 


If  the  instantaneous  frequency  is  a  linear  function  of  time  rather  than  some  arbitrary 
baseband  signal,  then  Eq.  (5.1)  is  referred  to  as  a  LFM  waveform.  Let  the  angle 
modulating  signal  be  of  the  form 

e(t)  =  Dpt2  (5.12) 

where  Dp  is  referred  to  as  the  phase-deviation  constant  or  phase  sensitivity.    Substituting 
Eq.  (5.12)  into  Eq.  (5.11)  results  in  the  following  expression  for  the  instantaneous 

frequency  : 


f i  (t)  =  f c  +  %^  .  (5.13) 


The  value  of  the  instantaneous  frequency  given  in  Eq.  (5.13)  is  a  linear  function  of  time  as 
desired,  composed  of  a  constant  carrier  frequency  and  the  term  (Dpt  /  7t)  which  is  known  as 
the  frequency  deviation.  The  value  (Dp  /  n)  in  this  instance  is  sometimes  referred  to  as  the 
frequency  sensitivity  of  the  modulator.  [Ref.  3:pp.  193-200]  [Ref.  12:pp.  267-271] 

Upon  substituting  Eq.  (5.12)  into  Eqs.  (5.1)  and  (5.6),  we  obtain  the  following 
expressions  for  the  LFM  waveform  and  its  complex  envelope  : 


xlfm  (t)  =  a(t)cos(27tfct  +  Dpt2)  (5.14) 


and 


2 
x1&n(t)  =  a(t)ejDPt    .  (5.15) 


On  comparing  the  above  expressions  with  Eqs.  (5.9)  and  (5.10),  it  can  be  seen  that  the 
expressions  for  the  CW  waveform  can  be  easily  obtained  from  those  for  the  LFM 
waveform  by  setting  the  phase  deviation  constant  equal  to  zero.  This  has  important 
consequences  for  the  computer  code  that  generates  these  waveforms.  It  means  that  we 
only  require  code  for  the  LFM  waveform  since  we  have  seen  that  CW  is  merely  a  special 
case  of  LFM. 

For  the  purposes  of  this  thesis  we  will  be  concerned  primarily  with  the  rectangular 
amplitude  modulating  function  given  by 

■»-tt:S>5 

where  A  is  a  constant  and  Tp  is  the  pulse  length  (in  seconds).  It  should  be  noted  from  Eq. 
(5.16)  that  the  pulse  is  assumed  to  be  centered  at  t  =  0,  and  that  the  modulating  function  is 
in  essence  acting  like  a  weighted  on/off  switch,  either  scaling  the  waveform  by  a  constant 
or  setting  it  equal  to  zero.  The  reason  that  we  are  interested  in  this  particular  amplitude 
modulating  function  is  that  it  is  because  it  corresponds  to  that  transmitted  by  a  "typical" 
SONAR  system.  Other  amplitude  modulating  functions  that  will  be  considered  are  the 
Hanning  and  Hamming  weighting  functions,  which  will  be  discussed  in  the  final  section  of 
this  chapter. 

For  a  LFM  pulse  with  a  pulse  length  of  Tp  (in  seconds),  the  quantity  (DpTp/7c)  is 
referred  to  as  the  swept  bandwidth  (SBW)  of  the  signal  x(t).  This  value  represents  the 
frequency  change  (in  Hz)  exhibited  by  the  signal  over  the  period  of  a  single  pulse  and 
assumes  that  only  one  sweep  occurs  during  the  period  of  the  pulse.  We  are  further 
assuming  that  this  swept  bandwidth  is  centered  around  the  carrier  frequency. 
[Ref.  3:pp.  197] 
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C.    FOURIER  SERIES  REPRESENTATION 

As  we  are  modelling  the  ocean  medium  in  terms  of  transfer  functions,  it  is  easier  and 
more  time  efficient  to  work  in  the  frequency  domain  as  opposed  to  the  time  domain. 
Accordingly,  we  require  x(t)  to  be  in  a  form  for  which  the  transform  X(f)  can  be  easily 
computed.  Except  for  a  few  simple  cases,  there  is  no  closed  form  expression  for  X(f) .  In 
fact,  X(f)  for  the  rectangular  envelope  LFM  pulse  is  similar  in  form  to  a  Fresnel  diffraction 
integral. 

Assume  that  we  can  expand  our  complex  envelope  into  the  truncated  complex  Fourier 
series 


i2rcqf0t        .  ..To 


x(t)=     I    cqej2TO,t°t       ;|t|<^  (5.17) 

q=-K 


where 


_TJ-To/ 


x(t)e-j27iqfol  dt  (5.18) 


is  the  complex  Fourier  series  coefficient  at  harmonic  q,  f0  =  (1/T0),  and  TD  is  sometimes 
referred  to  as  the  data  record  length  and  is  the  time  interval  over  which  we  are  attempting  to 
approximate  the  signal.  Since  this  series  is  expanded  about  a  period  of  TG,  the 
fundamental  frequency  of  the  Fourier  series  expansion  is  f0.  Accordingly,  the  coefficients 
represent  the  magnitude  and  phase  of  the  frequency  components  at  multiples  of  f0  Hz. 
This  can  be  clearly  demonstrated  by  taking  the  Fourier  transform  of  Eq.  (5.17)  which 
yields  the  result: 
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K 
X(f)=    X   cq8(f-qf0).  (5.19) 

q=-K 


The  value  K  determines  the  number  of  terms  contained  in  the  truncated  Fourier  series  and 
can  be  interpreted  as  the  number  of  harmonics  required  to  represent  the  complex  envelope 
x(t). 

The  effectiveness  of  the  series  representation  given  by  Eq.  (5.17)  lies  in  the  ability  of  a 
small  number  of  terms  to  provide  an  adequate  approximation  of  the  original  signal.  The 
complex  exponentials 


Si2",f°t  (5.20) 


used  in  this  series  expansion  form  an  orthogonal  basis  set  of  functions  [Ref.  15:pp.  214- 
215].  Consequently,  the  approximation  from  Eq.  (5.17)  using  the  coefficients  calculated 
from  Eq.  (5. 18)  results  in  the  minimum  mean  squared  error  estimate  of  x(t)  for  the  number 
of  terms  used  [Ref.  1 1  :pp.  32-36].  The  choice  of  complex  exponentials  as  the  basis  set  of 
functions  produces  a  series  expansion  that  relates  the  time  and  frequency  domains.  We 
have  to  consider  a  truncated  series  because  we  can  only  process  a  finite  amount  of  data,  and 
the  less  terms  we  need  to  represent  x(t)  the  less  computer  processing  will  be  required. 

The  preceding  discussion  of  the  Fourier  series  assumed  that  x(t)  was  a  continuous  time 
signal  but  the  concepts  and  results  can  equally  be  applied  to  a  discrete  time  version.  We 
will  now  discretize  x(t)  and  introduce  the  Discrete  Fourier  Transform  (  DFT{}).  The 
forward  and  inverse  DFT  can  be  defined  as 
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and 


DFT{x(l)}  =    £    xO)e-j27tql/L  (5.21) 

1=-L' 


IDFT{X(q)}=I   X   X(q)ej2m,1/L,  (5.22) 

q=-K 


respectively,  where  IDFT{}  is  the  Inverse  Discrete  Fourier  Transform  and  L  =  (2L'  +  1) 
is  the  total  number  of  time  samples  taken  of  x(t)  during  the  time  period  T0  .  To  digitize 
x(t),  we  will  sample  it  at  time  intervals  of  (TG/L)  which  results  in  x(l)  =  x(lTG/L). 
Substituting  this  digitized  version  of  x(t)  into  Eqs.  (5.17)  through  (5.19)  and  comparing 
with  Eqs.  (5.21)  and  (5.22),  it  can  be  shown  that 

cq=lX(q)  =  iDFT{x(l)}  (5.23) 


and 


x(l)=LxIDFT{cq}.  (5.24) 


Equations  (5.23)  and  (5.24)  in  conjunction  with  Eqs  (5.7)  and  (5.8)  are  the  primary 
equations  that  we  will  use  in  moving  between  the  time  and  frequency  domains  of  both  the 
bandpass  signal  and  its  complex  envelope. 

It  should  be  noted  that  this  concept  of  representing  a  signal  by  a  truncated  Fourier  series 
is  completely  general  and  can  be  applied  to  any  signal.  No  assumptions  regarding  the 
form  of  the  complex  envelope,  x(t),  were  made  in  the  development  of  this  section. 
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Consequently,  if  the  complex  envelope  is  of  the  form  given  by  Eq.  (5.6),  no  assumptions 
have  been  made  regarding  the  amplitude  and  angle  modulation  terms  of  Eq.  (5.1). 

D.    IMPLEMENTATION 

Let  us  now  consider  the  specific  case  of  rectangular  amplitude  modulated  CW  and  LFM 
pulses.  We  require  a  method  of  determining  the  value  of  K,  so  that  the  Fourier  series 
representation  contains  an  adequate  number  of  harmonics  to  accurately  represent  x(t),  and 
a  value  of  L  such  that  the  sampling  interval  is  sufficient  to  satisfy  Nyquist's  criterion  or 
some  other  specified  sampling  rate. 

The  complex  envelope  of  a  CW  pulse  can  be  thought  of  in  the  time  domain  as  a 
rectangular  window  convolved  with  a  sequence  of  delta  functions  at  intervals  of  TG 
seconds.  It  can  easily  be  seen  that  as  a  consequence,  the  resulting  DFT  will  be  a  discrete 
representation  of  the  sine  waveform.  To  represent  this  signal  we  will  take  all  the 
frequency  components  (separated  by  1/T0  Hz)  occurring  within  three  sine  function  zero 
crossings  (each  separated  by  1/Tp  Hz)  from  the  origin.  All  the  frequency  components 
extending  past  this  point  will  be  less  than  10%  of  the  value  of  the  frequency  component  at 
the  origin,  so  we  will  consider  them  as  being  insignificant  and  ignore  them,  which  yields 
the  relationship 


KCw=T-2-.  (5.25) 


As  stated  earlier,  there  is  no  closed  form  expression  for  the  LFM  pulse  in  the  frequency 
domain.  However,  for  a  rectangular  envelope  LFM  pulse  Papoulis  [Ref.  12:pp.  267-271] 
shows  that  if  the  following  condition  is  satisfied 
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(2fc-SBW)»V^  ,  (5.26) 

*      *  p 

then  the  resulting  transform  can  be  approximated  as  having  a  constant  magnitude  within  a 
frequency  band  equal  to  the  swept  bandwidth,  centered  about  the  carrier  frequency  fc  . 
Outside  of  this  region  the  transform  decays  to  zero.  We  need  to  approximate  this  rate  of 
decay  so  as  to  decide  how  many  frequency  components  we  need  to  accurately  represent  the 
original  bandpass  signal  x(t). 

The  corresponding  complex  envelope  of  the  LFM  pulse  will  therefore  be  approximately 
constant  within  the  region  -SBW/2  <  f  <  SBW/2  and  decaying  to  zero  outside  of  this 
region.  Since  the  frequency  components  are  spaced  1/T0  Hz  apart,  the  value  of  K  required 
to  represent  this  constant  region  of  the  spectrum  is  given  by 


K-  TEIo.  (5.27) 


Now  consider  the  CW  pulse  from  the  point  of  view  as  being  a  special  case  of  the  LFM 
pulse  where  the  swept  bandwidth  equals  zero.  We  would  therefore  expect  that  in  the  limit 
as  the  swept  bandwidth  goes  to  zero  that  Kifm  =  Kcw  .  We  will  therefore  approximate 
this  decaying  region  of  the  spectrum  as  being  the  width  of  three  zero  crossings  of  the  sine 
function  resulting  from  just  the  rectangular  amplitude  function.  It  is  surmised  that 
frequency  components  outside  this  region  will  be  less  than  10%  of  the  frequency 
component  at  the  origin  and  can  therefore  be  disregarded.  Therefore,  adding  the  K  values 
from  Eqs.  (5.25)  and  (5.27)  yields 


^SBWTo^  (52g) 
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which  is  valid  for  both  CW  and  LFM  pulses.  It  should  be  noted  from  Eq.  (5.17)  that  K  is 
an  integer  value,  so  in  our  computer  implementation  of  Eq.  (5.28),  K  is  rounded  up  to  the 
next  highest  integer  so  as  not  to  increase  the  truncation  error. 

Like  the  value  for  K,  L  is  also  an  integer,  but  with  the  additional  constraint  that  it  is  an 
odd  number.  This  is  in  order  to  satisfy  the  assumed  form  of  Eq.  (5.21)  and  to  make  the 
simulation  output  conform  to  that  required  by  the  adaptive  beamforming  algorithm  that  will 
be  used  to  validate  some  of  our  results.  The  number  of  samples  required  can  easily  be 
computed  from  the  calculated  values  of  K  using 

L  =  SR  x  K  (5.29) 

where  SR  is  the  user  defined  value  of  the  sampling  rate  and  K  is  the  value  determined 
from  Eq.  (5.28).  In  order  to  satisfy  Nyquist's  criterion,  the  sampling  rate  must  be  greater 
than  two.  This  means  that  the  sampling  frequency  must  be  greater  than  twice  the  highest 
frequency  component  present  in  the  truncated  series  representation  of  the  transmitted 
electrical  signal.  In  the  computer  implementation  of  Eq.  (5.29),  L  is  rounded  to  the  next 
highest  odd  integer  value  in  order  to  satisfy  all  of  the  previously  mentioned  constraints. 

In  the  computer  code,  the  amplitude  and  angle  modulated  signal  of  Eq.  (5.1)  is  defined 
by  the  following  variables  whose  values  are  supplied  by  the  user: 

•  AMP  :  the  magnitude  of  the  rectangular  amplitude  modulating  signal. 

•  TP  :  the  pulse  length  of  the  transmitted  pulse  (sec). 

•  PRF  :  the  pulse  repetition  frequency  (Hz).   This  value  is  the  same  as  fG  and  is 

equal  to  the  inverse  of  the  data  record  length  T0  (sec). 

•  FC  :  the  carrier  frequency  (Hz). 

•  SBW  :  the  swept  bandwidth  (Hz).    To  describe  a  CW  pulse,  set  this 

input  parameter  to  zero. 
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•  CHIRP  :  indicates  whether  the  LFM  sweep  is  increasing  in  frequency  (up)  or 

decreasing  (down). 

•  SR  :  sampling  rate.    This  value  is  the  ratio  of  the  sampling  frequency  to  the 

maximum  frequency  component  present  in  the  signal  and  must  be  greater 
than  two  in  order  to  satisfy  Nyquist's  criterion. 

•  AMOD  :  the  amplitude  modulation  function,  selected  from  a  menu  of  choices.  The 

rectangular  amplitude  modulation  function  is  denoted  by  AMOD  =  0. 


E.    THE  RECTANGULAR  AMPLITUDE  MODULATED  PULSE 

To  validate  the  theory  of  the  preceding  section,  examples  of  rectangular  amplitude 
modulated  CW  and  LFM  pulses  were  reconstructed  from  the  calculated  Fourier  series 
coefficients  of  its  complex  envelope,  denoted  xce(n),  and  compared  with  the  actual 
sampled  bandpass  signal  denoted  x(n).  In  the  following  plots,  x(n)  is  plotted  as  a  dashed 
curve  while  Xce(n)  is  plotted  as  a  solid  curve.  For  plotting  purposes,  a  third-degree 
polynomial  curve  fitting  routine  is  used  to  interpolate  between  the  data  points. 
1.    CW  Pulse  with  50%  Duty  Cycle 

Duty  cycle  is  defined  as  the  ratio  Tp/TG,  which  for  this  case  means  that  the 
amplitude  weighting  function  is  nonzero  for  50%  of  the  time.  While  this  is  not  a 
particularly  practical  transmit  signal  for  underwater  acoustic  applications,  it  will 
demonstrate  that  a  real  bandpass  signal  can  be  reconstructed  from  its  complex  envelope. 
Combining  this  high  value  of  duty  cycle  with  a  suitably  low  carrier  frequency,  the 
comparative  plots  of  xce(n)  and  x(n)  will  be  uncluttered  and  will  help  illustrate  some 
important  points. 

a.  Pulse  Reconstruction 

A  CW  pulse  with  the  following  characteristics  is  shown  in  Fig.  5.1a 

•     AMP  =  1,  TP  =  0.0125  sec,  PRF  =  40  Hz,  FC  =  400  Hz,  SR  =  4.0 
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It  can  be  seen  from  Fig.  5.1a  that  the  two  curves,  x(n)  and  xce(n),  agree 
fairly  closely.  This  indicates  that  we  are  indeed  able  to  reconstruct  a  real  bandpass  signal 
from  its  complex  envelope.  However,  there  are  some  important  effects  to  be  observed. 
Note  the  manner  in  which  xce(n)  overshoots  x(n)  when  the  pulse  is  turned  on  and  off.  It 
can  also  be  seen  that  xce  (n)  tends  to  oscillate  around  the  true  value  of  x(n). 

A  truncated  Fourier  series,  by  its  very  nature,  will  exhibit  high  frequency 
oscillations  about  the  true  value  of  the  function.  While  these  oscillations  usually  have 
sufficiently  small  amplitudes  to  be  of  no  major  consequence,  this  is  not  the  case  at  a 
discontinuity.  Where  x(n)  is  turned  on/off  corresponds  to  just  such  a  discontinuity  in  the 
rectangular  amplitude  weighting  function.  As  more  terms  are  added  to  the  truncated  series 
the  overshoot  does  not  tend  to  disappear.  Instead,  it  becomes  narrower  due  to  the  higher 
frequency  terms  present  and  the  first  overshoot  approaches  a  value  of  0.08949  of  the 
discontinuity.  This  behavior  is  known  as  "Gibbs  Phenomenon"  and  is  due  to  the  truncated 
Fourier  series  representation  of  the  complex  envelope,  rather  than  the  fact  that  we  have 
used  the  complex  envelope  to  reconstruct  the  original  bandpass  signal. 
[Ref.  13:pp.  107-112]  [Ref.  14:pp.  179-185] 
b.    Lanczos  Smoothing 

It  was  suggested  by  Cornelius  Lanczos  [Ref.  15:pp.  219-221]  that  this 
phenomenon  could  be  reduced  if  the  Fourier  series  could  be  made  to  converge  more 
quickly.  He  showed  that  if  the  series  coefficients  were  modified  by  a  correction  term  he 
called  the  sigma  factors,  given  by 


sin^q7p/K) 
0(K'q)"    (q7t/K)    .  (530) 
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where  K  is  the  order  of  the  last  term  in  the  series,  then  the  tendency  of  the  high  order  terms 
to  make  a  series  divergent  could  be  counteracted.  As  can  be  seen  from  Eq.  (5.30)  this  was 
achieved  by  increasing  the  amount  of  attenuation  with  increasing  order  of  the  Fourier 
coefficients.  Hamming  [Ref.  13:pp.  107-112]  comments  that  K  in  Eq.  (5.30)  is  usually 
replaced  with  (K+l),  and  this  is  the  convention  we  will  adopt  in  our  application  of  the 
sigma  factors  for  this  thesis. 

Lanczos  also  showed  that  the  application  of  the  sigma  factors  was 
analogous  to  averaging  over  an  incremental  time  period  of  ±T0/  2K.  This  means  that  each 
original  point  in  the  truncated  Fourier  series,  x(t),  is  replaced  by  the  arithmetic  mean  around 
that  point,  between  the  limits  of  the  incremental  time  period.  Since  the  ripple  in  the 
truncated  series  representation  has  the  period  of  either  the  first  term  neglected  or  the  last 
term  kept,  the  main  effects  of  this  ripple  can  be  removed  by  integrating  or  averaging  over 
this  period.  The  averaging  interval  usually  chosen  for  this  smoothing  is  TG  /(K+l),  which 
corresponds  to  replacing  K  with  (K+l)  in  Eq.  (5.30).  If  K  is  large,  then  the  region  of 
local  averaging  is  very  small  and  has  very  little  effect,  except  in  the  neighborhood  of  the 
discontinuity.  This  is  a  desirable  result  since  the  discontinuity  is  the  region  where  the 
smoothing  is  most  needed. 

When  applied  to  the  truncated  Fourier  series,  this  is  known  as  "Lanczos 
Smoothing".  Applying  this  smoothing  to  the  truncated  series  of  Eq.  (5.17)  yields  the 
following: 


s(t)=      £    o(K,q)c  c 


.J27cqf0t 


;M* 


(5.31) 


As  can  be  seen  from  Eq.  (5.31),  the  smoothed  series  is  nothing  more  than  the  original 
Fourier  series  with  its  coefficients  multiplied  by  the  corresponding  sigma  factors.  This 
result  is  what  we  would  expect  owing  to  the  linearity  of  the  Fourier  series.  The  effect  of 
time  averaging  is  the  same  as  convolving  the  time  sequence  with  a  weighted  rectangular 
window.  This  in  turn  is  analogous  to  multiplying  the  corresponding  Fourier  series 
coefficients  of  the  time  sequence  with  those  of  the  window.  Thus,  the  sigma  factors  can 
be  thought  of  as  the  Fourier  series  coefficients  of  this  weighted  rectangular  window. 
c.    The  Smoothed  Pulse 

The  effect  of  the  application  of  these  sigma  factors  can  be  seen  in  Fig.  5.1b, 
where  the  Lanczos  Smoothing  has  been  applied  to  the  Fourier  Series  coefficients  of  our 
CW  pulse.  As  expected,  the  oscillations  and  overshoot  at  the  discontinuity  have  been 
reduced  but  at  the  expense  of  slightly  broadening  the  pulse  and  reducing  the  rate  of  change 
at  the  leading  and  trailing  edges  of  the  pulse.  It  should  be  noticed  that  the  oscillations  have 
been  reduced  for  both  portions  of  the  duty  cycle,  i.e.,  not  just  during  the  pulse. 

The  first  effect  is  to  be  expected  since  we  are  effectively  convolving  our 
time  sequence  with  a  windowing  function.  Whenever  we  convolve  two  sequences  the 
result  is  a  sequence  larger  than,  but  related  to,  the  size  of  the  original  functions.  Since  the 
windowing  function  is  significantly  smaller  than  the  pulse  duration,  we  expect  the  effects 
of  this  spreading  to  be  minimal.  On  comparing  Figs.  5.1a  and  5.1b  we  can  see  that  this  is 
certainly  the  case  for  this  example. 

The  second  effect  is  due  to  the  fact  that  the  smoothing  technique  attenuates 
the  high  order  Fourier  coefficients.  Attenuating  high  order  coefficients  is  analogous  to 
attenuating  the  high  frequency  components  of  the  series  representation  of  the  bandpass 
signal.  Thus,  we  would  expect  the  smoothed  series  representation  of  the  signal,  xce(n),  to 
be  slower  than  the  original  Fourier  series  at  reacting  to  sudden  change,  i.e.,  a  discontinuity. 
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-12.5     -10.0       -7.5       -5.0       -2.5        0.0        2.5 

TIME  T(MSEC) 
CRSE: FIG  5.1H  WRVEF0RM:CW  PULSE  I  NFREQ:  13  ) 
R:  1.0   TF:  12.5  MSEC   FRF:  40.0  HZ   FC=  400.0  HZ 


Fig.  5.1a     CVV  pulse,  50%  duty  cycle,  no  smoothing. 


-12.5     -10 


TIME  T(MSEC) 

CASE: FIG  5. IB     WAVEFORM :CW  TUL5E  (  NrREQ:  13  ) 
R:  1.0   TP:  12.5  MSEC   FRF:  40.0  HZ   FC:  400.0  HZ 


0.0  12.5 


Fig.  5.1b     CVV  pulse,  50%  duty  cycle,  Lanczos  smoothing. 
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Regarding  the  application  of  the  sigma  factors  to  the  specific  problem  of  Gibbs 
phenomenon,  Lanczos  [Ref.  15:pp.  225-226]  notes  that  unlike  the  often-used  Fejer 
arithmetic  mean  method,  the  sigma  factors  do  not  completely  eliminate  the  oscillations  due 
to  the  discontinuity.  However,  the  oscillations  are  significantly  reduced  and  the  sigma 
factors  do  not  cause  nearly  as  drastic  a  reduction  in  the  rate  of  change  at  the  leading  and 
trailing  edges  of  the  pulse.  It  should  be  noted  that  this  smoothing  invalidates  the 
minimum-mean-square-error  property  of  the  Fourier  series.  The  original  Fourier  series 
coefficients  are  the  only  ones  that  satisfy  this  property  and  now  they  have  been  changed. 

2.    LFM  Pulse  with  50%  Duty  Cycle 

A  LFM  pulse  with  the  following  characteristics  is  shown  in  Fig.  5.2a 


AMP  =  1,  TP  =  0.0125  sec,  PRF  =  40  Hz,  FC  =  400  Hz,  SRATE  =  4.0 
SBW  =  60Hz,  "up"  chirp 


The  purpose  of  Figs  5.2a  through  5.2c  is  to  demonstrate  that  our  analysis  of  the  LFM  pulse 
was  correct.  As  can  be  seen  from  Fig.  5.2a,  we  are  able  to  reconstruct  the  LFM 
waveform  from  the  Fourier  series  representation  of  its  complex  envelope,  as  was  the  case 
for  the  CW  pulse.  The  LFM  pulse  in  this  case  is  said  to  have  "up  chirp"  because  its  phase 
deviation  constant  Dp  is  positive,  that  is,  the  frequency  is  increasing  during  the  period  of 
the  pulse. 

As  for  the  CW  pulse,  the  reconstructed  signal  xce(n)  oscillates  about  its  true 
value  x(n)  and  exhibits  the  characteristic  overshoots  at  the  leading  and  trailing  edges  of  the 
pulse.  The  application  of  the  sigma  factors  in  Fig.  5.2b  shows  that  the  Lanczos  smoothing 
is  as  effective  for  the  LFM  case  as  it  was  for  the  CW  pulse  and  exhibits  the  same  results. 


64 


TIME  TIMSEC) 
CP5E:  FIG  5.2A     HRVEFORH:LFH  PULSE  (  NFREO:  15  ) 
Rs  1.0   TP:  12.5  MSEC   FRF:  40.0  HZ   FC  =  100.0  HZ   SWPTBW:  60.0  HZ   CHIRP:  U 

Fig.  5.2a     LFM  pulse,  50%  duty  cycle,  "up  chirp",  no  smoothing. 


TIME  TIMSEC) 
CASE: FIG  5.25     WRVEF ORH: LFM  PULSE (  NTREOs  15  ) 

R:  1.0   TP:  12.5  MSEC   PRF:  40.0  HZ   FC:  400.0  HZ   SWPTBW:  60.0  HZ   CHIRP:  U 

Fig.  5.2b     LFM  pulse,  50%  duty  cycle,  "up  chirp",  Lanczos  smoothing. 
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It  sh6uld  be  noted  however,  that  since  the  swept  bandwidth  is  low  and  the  pulse 
repetition  frequency  is  fairly  high,  the  value  of  K  given  by  Eq.(5.27)  is  less  than  one. 
This  effectively  means  that  the  signal  is  approaching  the  limit  of  becoming  a  CW  pulse. 
Inspection  of  Eq.  (5.28)  shows  that  the  dominating  effect  on  the  number  of  significant 
frequency  components  is  due  to  the  rectangular  pulsed  nature  of  the  waveform  rather  than 
the  frequency  modulation.  Further  investigation  needs  to  be  carried  out  on  whether  or  not 
Eq.  (5.27)  is  a  valid  approximation.  This  will  be  done  by  the  study  of  a  LFM  pulse  with  a 
smaller  duty  cycle  to  follow. 

The  LFM  pulse  in  Fig.  5.2c  differs  from  the  previous  LFM  pulses  in  that  it  is  a 
"down  chirp"  pulse,  that  is  its  phase  deviation  constant  is  negative  and  the  frequency  is 
decreasing  during  the  period  of  the  pulse.  The  purpose  of  this  figure  is  to  show  that  the 
computer  code  can  model  the  "down  chirp"  pulse  as  well  as  the  "up  chirp"  case. 


TIME  T(HSEC) 
CRSE:FIG  5.2C     WRVEFORH : LFM  PULSE (  NFREO: 15  ) 
R:  1.0   TP:  12.5  MSEC  PRF:  40.0  HZ  FC:  400.0  HZ  SWPIBW:  60.0  HZ  CHIRP:  D 

Fig.  5.2c     LFM  pulse,  50%  duty  cycle,  "down  chirp",  Lanezos  smoothing. 
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3.    CW  PULSE  WITH  8%  DUTY  CYCLE 

A  pulse  waveform  with  a  small  duty  cycle  is  of  more  practical  use  for 
underwater  acoustics  applications,  so  we  need  to  check  that  our  computer  code  will  handle 
more  realistic  signals.  A  CW  pulse  with  the  following  characteristics  is  shown  in  Figs. 
5.3a  and  5.3b: 


•     AMP  =  1,  TP  =  0.04  sec,  PRF  =  2.0  Hz,  FC  =  400  Hz,  SR  =  4.0 

As  can  be  seen  from  the  above  parameters,  the  ratio  Tp/T0  yields  the  value  0.08  which 
corresponds  to  a  duty  cycle  of  8%.  In  Figs.  5.3a  and  5.3b,  only  the  curve  reconstructed 
from  its  Fourier  series  coefficients,  xce(n),  is  shown  so  that  the  plots  do  not  appear  too 

cluttered. 

From  Fig.  5.3a  we  can  observe  that  as  expected,  the  unmodified  Fourier 
coefficients  produce  a  waveform  that  oscillates  about  the  true  value  of  x(n).  The  ringing, 
or  oscillation,  at  the  leading  and  trailing  edges  of  the  pulse  is  quite  pronounced  and  give  the 
appearance  of  the  rectangular  pulse  shape  spreading. 

The  application  of  Lanczos  smoothing  to  this  case,  shown  in  Fig.  5.3b, 
demonstrates  that  the  oscillations  about  the  true  value  have  been  reduced  but  not  completely 
removed.  This  can  be  seen  to  be  at  the  expense  of  a  slight  spreading  of  the  pulse  and  a 
decrease  in  the  rate  of  change  at  the  leading  and  trailing  edges  of  the  pulse.  As  expected, 
these  are  the  same  effects  as  were  observed  for  the  CW  pulse  with  the  50%  duty  cycle 
except  that  the  extent  of  the  ringing  present  in  this  case  is  more  dramatic. 
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TIME  T(MSEC) 
CflSE:FIG  5.3R  WAVEFORM :CW  PULSE  (  NFREO: 77  ) 
Hi  1.0   TP:  -10.0  MSEC  PRF:  2.0  HZ   FC:  400.0  HZ 


Fig.  5.3a     CW  pulse,  8%  duty  cycle,  no  smoothing. 


-250.0        -200.0  -150. 0  -100. 0     .        -50.0  0.0 

TIME  T(MSEC) 
CflSEsFIQ  5.3B  WRVEFORH:CW  PULSE t   NFREO: 77    ) 

H:    1.0      TP:    40.0    MSEC      PRF:    2.0    HZ     FC:    400.0    HZ 


Fig.  5.3b     CW  pulse,  8%  duty  cycle,  Lanczos  smoothing. 
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Why  are  we  interested  in  this  amount  of  ringing  in  the  pulse  representation? 
One  of  the  outputs  from  this  computer  simulation  will  be  the  time-domain  output  electrical 
signal  at  each  element  in  the  receive  array.  We  wish  to  use  this  information  to  illustrate  the 
pulse  distortion  due  to  dispersion  in  the  ocean  waveguide.  Consequently,  we  wish  to 
"clean  up"  the  transmitted  pulse  as  much  as  possible  so  that  we  do  not  confuse  distortion 
due  to  the  Fourier  series  representation  of  the  transmitted  signal  with  distortion  due  to 
dispersion  in  the  ocean  waveguide. 

By  comparing  Figs.  5.3a  and  5.3b,  it  can  be  seen  that  the  application  of 
Lanczos  smoothing  does  produce  a  much  "cleaner"  pulse  shape.  It  is  our  contention  that  it 
will  be  easier  to  observe  the  effects  of  dispersion  on  this  smoothed  pulse  that  on  the 
unmodified  pulse. 

4.    LFM  Pulse  With  8%  Duty  Cycle 

A  LFM  pulse  with  the  following  characteristics  is  shown  in  Figs.  5.4a  and 
5.4b: 


AMP  =  1,  TP  =  0.04  sec,  PRF  =  2.0  Hz,  FC  =  400  Hz,  SRATE  =  4.0 

SBW  =  120  Hz,  "up  chirp" 


As  for  the  CW  case,  we  have  only  presented  the  signal  reconstructed  from  its  Fourier  series 
coefficients,  Xce(n),  so  that  the  plots  do  not  become  too  cluttered. 

It  can  be  observed  that  all  the  comments  made  regarding  the  CW  pulse  with  the 
8%  duty  cycle  are  applicable  here.  It  can  be  seen  from  comparing  Figs.  5.3a  and  5.4a 
however,  that  the  oscillations  about  the  true  value  and  the  amount  of  ringing  appears  to  be 
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-100.0  -50.0  0.0  50.0  100.0  150.0  200.0 

TIME  TU1SEC) 
CnSE:FIG5.4fl  WflVEFORH :  LFM  PULSE  (    NFREO:  137    ) 

fi:     1.0      TP:    40.0    MSEC      FRF:    2.0    HZ      FC:    -500.0    HZ      SWFTBW:     120.0    HZ      CHIRP:? 

Fig.  5.4a     LFM  pulse,  8%  duty  cycle,  "up  chirp",  no  smoothing. 


-200.0  -150.0  -1D0.0  -50.0  0  0  50.0  100J0  150.0  200.0 

TIME   T(MSEC) 

:  FIG  5.4B  WFlVEFORM :  LFM  PULSE  (    NFREO:  137  ) 

.0      TP:    40.0    MSEC      FRF:    2.0    HZ      FC:    400.0  HZ      SWP1BW:     120.0    HZ      CHIRP:  7 


Fig.  5.4b     LFM  pulse,  8%  duty  cycle,  "up  chirp",  Lanczos  smoothing. 
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less  for  the  LFM  pulse  than  it  does  for  the  CW  pulse.    This  is  probably  due  to  the  number 
of  frequency  components  used  to  represent  the  signal. 

It  would  appear  that  the  value  of  K  given  by  Eq.  (5.28)  is  a  more  conservative 
estimate  of  the  number  of  frequency  components  required  to  represent  the  LFM  pulse  than 
the  estimation  for  the  CW  pulse  given  by  Eq.  (5.25).  As  the  number  of  frequency 
components  used  is  increased,  the  amount  of  ringing  at  the  leading  and  trailing  edges 
should  be  reduced.  Since  the  amount  of  ringing  for  the  LFM  waveform  is  less  than  for  the 
CW  case,  we  assume  that  we  erred  on  the  side  of  caution  and  included  more  frequency 
components  that  absolutely  necessary  to  represent  a  recognizable  LFM  pulse  shape.    This 

is  consistent  with  the  logic  we  used  in  proposing  the  form  of  Eq.  (5.28) 

This  would  also  explain  why  the  application  of  the  Lanczos  smoothing  in  Fig. 
5.4b  appears  to  have  distorted  the  LFM  pulse  shape  slightly  less  than  it  did  for  the  CW 
pulse  in  Fig.  5.3b.  If  we  are  at  the  minimum  number  of  frequency  components  required  to 
represent  the  signal,  then  the  attenuation  of  the  high  order  terms  will  have  a  more  dramatic 
effect  than  if  we  have  a  slight  excess  in  the  number  of  frequency  components  required. 
This  effect  can  be  shown  quite  dramatically  in  Figs.  5.4c  and  5.4d. 

In  Figs.  5.4c  and  5.4d,  we  have  reconstructed  the  LFM  pulse  using  the  number 
of  frequency  components  given  by  Eq.  (5.27),  which  assumes  that  only  the  constant 
portion  of  the  spectrum  of  the  LFM  pulse  is  significant.  As  can  be  seen  from  Fig.  5.4c, 
the  reconstructed  pulse  xce(n)  is  far  from  being  a  clean  looking  pulse  shape.  It  oscillates 
erratically  about  the  true  value  and  the  amount  of  ringing  exhibited  is  quite  extensive. 
When  we  apply  the  Lanczos  smoothing  to  this  case  the  waveform  appears  distorted  as 
expected.  In  fact,  it  is  distorted  to  the  extent  that  we  are  not  able  to  tell  by  looking  at  it  that 
it  has  a  rectangular  pulse  shape. 


if 


-250.0        -200. 0  -150.0  -100.0  -50.0  00  50.0  100.0  150. 0  200.0 

TIME  T(NSEC) 
CR5E:FIG5.4C  UnVErCRM : LFM  PULSE (    MFRLQ: 63    ) 

fl:     1.0      TP:    40.0    MSEC      FRF':     2.0    HZ      FC:    -500.0    HZ      SWPTBH:     120.0    HZ      CHIRP:  U 

Fig.  5.4c     LFM  pulse,  8%  duty  cycle,  "up  chirp",  no  smoothing, 
insufficient  frequency  components. 
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200.0 


-250.0        -200.0  -150.0  -100.0  -50.0  0.0  50.0 

TIME  T(MSFC) 
Cn5E:FIG5.4D  WnVEFORM:LFM  FULSE  (    NTREQ: 63    ) 

fl:     1.0      TP:    -50.0    MSEC      FRF:    2.0    HZ      FC:    400.0    HZ      5WF1BW:     120.0    HZ      CHIRP: U 

Fig.  5.4d     LFM  pulse,  8%  duty  cycle,  "up  chirp",  Lanczos  smoothing, 
insufficient  frequency  components. 
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Clearly,  the  assumption  that  only  the  constant  portion  of  the  spectrum  of  the 
LFM  pulse  is  significant  is  not  valid.  Inspection  of  Figs.  5.4a  through  5.4d  justifies  our 
scheme  to  calculate  the  number  of  frequency  components  required  to  represent  a  LFM  pulse 
shape  given  by  Eq.  (5.28).  Note  that  Eq.  (5.28)  was  derived  for  the  case  of  rectangular 
amplitude  modulated  pulses.  Since  the  rectangular-shaped  function  contains  a  severe 
discontinuity,  we  expect  that  the  number  of  frequency  components  given  by  Eq.  (5.28)  will 
coincide  to  a  worst  case  situation.  Consequently,  we  would  expect  that  if  we  chose  a 
modulating  function  that  did  not  have  as  large  a  discontinuity  at  the  leading  and  trailing 
edges  of  the  pulse,  we  would  not  require  as  many  frequency  components  to  accurately 
represent  the  signal. 

F.   ALTERNATIVE  AMPLITUDE  MODULATION  FUNCTIONS 

As  discussed  earlier  in  this  chapter,  although  we  are  primarily  interested  in  the 
rectangular  amplitude  modulating  function,  it  is  not  the  only  one  that  will  be  considered. 
Of  interest  is  the  Hanning  weighting  function  given  by 


rA[0.5  +  0.5cos(27n/Tp)]  ,  lt|<Tp/2  (5  32) 


and  the  Hamming  weighting  function  given  by 


r  A[0.54+0.46cos(27rt/Tp)]  ,  |t|<Tp/2  (5  33) 


where  A  is  a  constant  and  Tp  is  the  pulse  length  (in  seconds).  It  should  be  noted  from 
Eqs.  (5.32)  and  (5.33)  that  the  pulse  is  assumed  to  be  centered  at  t  =  0,  as  was  the  case  for 
the  rectangular  amplitude  modulating  function  given  by  Eq.  (5.16).    Note  that  Eqs.  (5.32) 
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and  (5.33)  are  the  actual  pulse  amplitude  modulation  functions  and  should  not  be  confused 
with  the  concept  of  windowing  functions  applied  to  a  block  of  data,  i.e.,  the  modulation 
functions  are  applied  specifically  to  the  pulse  (of  length  Tp  seconds)  and  not  to  the  entire 
data  record  (of  length  TQ  seconds). 

Inspection  of  Eqs.  (5.32)  and  (5.33)  reveals  that  they  are  both  a  maximum  at  the  center 
of  the  pulse,  and  decrease  towards  zero  at  the  leading  and  trailing  edge  of  the  pulse  in  a 
"smooth"  fashion.  It  can  also  be  readily  seen  that  these  two  equations  are  very  nearly 
identical.  Why  then  are  we  interested  in  two  very  similar  amplitude  modulation  functions 
of  this  form?  From  the  discussion  in  the  previous  section,  we  expect  that  it  will  require 
less  frequency  components  to  accurately  represent  the  transmitted  signal  than  was  the  case 
when  a  rectangular  amplitude  modulation  function  was  used.  To  test  this  hypothesis, 
examples  of  amplitude  modulated  C W  pulses  are  plotted  as  follows: 

•  Fig.  5.5a.  Hanning  amplitude  modulation,  50%  duty  cycle, 

•  Fig.  5.5b.  Hanning  amplitude  modulation,  8%  duty  cycle, 

•  Fig.  5.6a.  Hamming  amplitude  modulation,  50%  duty  cycle, 

•  Fig.  5.6b.  Hamming  amplitude  modulation,  8%  duty  cycle. 

In  Figs.  5.5a  and  5.6a,  the  pulse  reconstructed  from  the  calculated  Fourier  series 
coefficients  of  its  complex  envelope,  denoted  xce(n),  is  compared  with  the  actual  sampled 
bandpass  signal,  denoted  x(n),  where  x(n)  is  plotted  as  a  dashed  curve  while  xce(n)  is 
plotted  as  a  solid  curve.  For  plotting  purposes,  a  third-degree  polynomial  curve  fitting 
routine  is  used  to  interpolate  between  the  data  points.  In  Figs.  5.5b  and  5.6b,  just  Xce  (n) 
is  plotted  to  prevent  the  plot  from  becoming  too  cluttered. 
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Comparing  Figs.  5.5a  and  5.6a  with  Fig.  5.1a,  it  can  be  seen  that  indeed  fewer 
frequency  components  are  required  to  accurately  represent  the  transmitted  signal  when  the 
amplitude  modulation  function  is  of  the  form  given  by  Eqs.  (5.32)  or  (5.33),  than  if  it  is  of 
the  form  given  by  Eq.  (5.16).  Whereas  we  introduced  the  concept  of  Lanczos  smoothing 
for  the  rectangular  modulated  pulse  in  an  attempt  to  obtain  a  "cleaner"  representation  of  the 
signal,  we  find  that  this  is  totally  unnecessary  for  the  Hanning  and  Hamming  modulation 
functions.  The  observations  that  can  be  made  when  comparing  Figs.  5.5b  and  5.6b  with 
Fig.  5.3a  are  consistent  with  those  just  described.  From  these  comparisons  we  can  make 
the  statement  that  the  Hanning  and  Hamming  modulated  CW  pulses  require  roughly  half 
the  number  of  frequency  components  for  accurate  representation  than  does  the  rectangular 
modulated  pulse.  No  comparison  of  modulated  LFM  pulses  was  carried  out  as  it  was 
judged  past  the  scope  of  interest  for  this  thesis. 

Further  comparison  of  Figs.  5.5  and  5.6  shows  that  for  a  given  number  of  frequency 
components  the  representation  of  the  Hamming  amplitude  modulated  pulse  is  better  than  the 
Hanning  modulated  pulse.  This  means  that  there  is  closer  agreement  between  the  pulse 
reconstructed  from  the  calculated  Fourier  series  coefficients  of  its  complex  envelope  and  the 
actual  sampled  bandpass  signal  when  the  Hamming  amplitude  modulation  function  is  used. 
In  particular,  the  Hamming  modulated  pulse  displays  significantly  less  ringing  before/after 
the  leading/trailing  edges  of  the  pulse  for  the  lower  duty  cycle  case  shown  in  Figs.  5.5b 
and  5.6b. 

The  advantage  in  reducing  the  number  of  frequency  components  required  to  represent  a 
signal  is  a  simple  one.  Since  the  computer  simulation  must  compute  a  value  of  the  overall 
system  complex  frequency  response  for  each  frequency  of  interest,  and  this  operation 
comprises  the  majority  of  the  computational  effort,  a  decrease  in  the  number  of  required 
frequencies  will  mean  an  almost  proportional  drop  in  computation  time,  i.e.,  if  the 
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TIME  T!MSEC> 
CnSE:FIG5.5n  WflVEFORHsCW  PULSE  (    NFREQsS    ) 

fl:     1.0      TF:     12.5    MSEC      FRF:     40.0    HZ      FC:     100.0    HZ 

Fig.  5.5a     CW  pulse,  50%  duty  cycle,  Manning  amplitude  modulation. 
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TIME  T(MSEC) 
CF1SE:  FIG  5.5B  WFV.'EFORIhCW  PULSE  (    NFREQ:  39    ) 

fl:     1.0      TP:    10.0    MSEC      FRF:    2.0    HZ      FC:    400.0    HZ 
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Fig.  5.5b     CW  pulse,  8%  duty  cycle,  Hanning  amplitude  modulation. 
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FC:    100.0    HZ 

Fig.  5.6a     C\V  pulse,  509c  duty  cycle,  Hamming  amplitude  modulation. 
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Fig.  5.6b     C\V  pulse,  8%  duty  cycle,  Hamming  amplitude  modulation. 
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number  of  frequencies  is  reduced  by  half,  then  the  program  run  time  will  be  approximately 
reduced  by  half. 

However,  this  is  not  the  main  reason  that  we  consider  amplitude  modulation  functions 
of  this  form.  As  stated  earlier,  we  are  primarily  interested  in  the  rectangular  amplitude 
modulating  function.  Our  interest  in  the  Hamming  and  Hanning  functions  stems  from  the 
fact  that  many  currently  available  computer  models  dealing  with  sound  propagation  in  the 
ocean  deal  with  transmitted  waveforms  having  these  particular  amplitude  modulation 
functions,  presumably  in  an  attempt  to  reduce  computation  times.  If  we  then  also  have  the 
option  of  using  these  same  modulation  functions,  we  will  be  in  a  better  position  to  compare 
our  results  and  performance  with  other  ocean  medium  propagation  models,  e.g.,  SAFARI 
[Ref.  7]  . 


78 


VI.   COMPUTER  SIMULATION  RESULTS  FOR  THE 
ISOSPEED  MEDIUM  CASE 

A.  THE  OVERALL  SYSTEM  COMPLEX  FREQUENCY  RESPONSE  FOR  AN 
ISOSPEED  MEDIUM 

The  first  step  in  validating  the  developments  carried  out  in  Chapters  2  through  4  is  to 
apply  the  analysis  to  the  simplest  scenario,  i.e.,  the  free-space,  isospeed  medium  case. 
Since  the  preceding  development  of  the  overall  system  complex  frequency  response  given 
by  Eq.  (3.14)  assumed  a  free-space  problem,  all  we  need  to  do  is  apply  the  condition  of  an 
isospeed  medium  to  the  ocean  medium  transfer  function  given  by  Eq.  (4.34).  As 
previously  discussed,  the  phase  integral  of  the  ocean  medium  transfer  function  has  a  closed 
form  expression  in  the  isospeed  case,  which  results  in  the  relationship 

HM(f.fr.yo;y)=  ^M^e^W,.)/^^*,,, 

2kY(y0)V|kY(y)|  v     ' 


where 


2      2    1/2        2  2 

kY(y0)  =  ±2*[(f/c0)  -rT]      ,    fr<(f/c0),  (6.2a) 


2  l}11        2  2 

kY(y0)  =  +  j27t[f; r-(f/c0n      ,    fr>(f/c0),  (6.2b) 


kY(y0)  =  27cfY(y0)  =  27cfY,  and  c0  is  the  constant  speed  of  sound  in  the  isospeed 
medium.     Equation  (6.1),  then,  is  our  ocean  medium  transfer  function  for  a  free-space, 
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isospeed  medium.  Note  that  the  ocean  medium  transfer  function  given  by  Eq.  (6.1) 
contains  no  integral  operations  due  to  there  being  a  closed  form  expression  for  the  phase 
integral  of  Eq.  (4.34).  This  result  is  significant  since  numerical  integration  routines  are 
very  time  consuming.  Therefore,  we  expect  that  the  computational  effort  required  for  the 
isospeed  case  will  form  a  baseline  against  which  to  compare  other  cases  that  do  not  have 
closed  form  solutions. 

It  can  be  easily  seen  from  Eqs.  (6.2a)  and  (6.2b)  that  for  an  isospeed  medium,  the 
propagation  vectors  in  the  Y  directions,  k Y  (yn )  and  k Y  (y),  are  equal.  The  ray  acoustics' 
interpretation  of  this  is  that  rays  travel  in  straight  lines  and  will  not  bend.  It  should  be 
remembered  from  the  development  in  Chapter  3,  that  Eq.  (6.2a)  represents  propagating 
waves  while  Eq.  (6.2b)  represents  evanescent  waves. 

Returning  to  our  expression  for  the  overall  system  complex  frequency  response  given 
by 

r°°      MT'  NT 

H(f,r)  =  27c   I  X  X       cmn(OHM(f,fr,yo;y)Jo(27tfrrm) 

•'O      m=-MT    n  =  -NT 

x  eJ^Y(y-y„)frdfr  (314) 

we  notice  that  the  region  of  integration  with  respect  to  fr  is  from  zero  to  infinity.  Equation 
(3.14)  is  referred  to  as  a.  full-wave  solution  for  our  ocean  medium  propagation  problem, 
and  can  be  thought  of  as  an  infinite  summation  of  ray  acoustic  solutions  [Ref.  7:p.  36]. 
From  Eqs.  (6.2a)  and  (6.2b),  it  can  be  seen  that  this  region  includes  both  propagating  and 
evanescent  waves.  Since  the  propagation  vector  in  the  Y  direction  has  different  forms 
depending  on  whether  the  wave  is  propagating  or  not,  we  will  split  Eq.(4.33)  into  two 
integrals.     One  integral  will  cover  the  region  of  propagating  waves,  i.e.,  integrate  with 


respect  to  fr  from  zero  to  (f  /  cG),  while  the  second  will  cover  the  region  of  evanescent 
waves,  i.e.,  integrate  with  respect  to  fr  from  (f  /  cQ)  to  infinity. 

As  discussed  in  Chapter  4,  our  general  solution  for  the  case  of  a  unit  amplitude  point 
source  in  an  unbounded,  homogeneous  medium  should  collapse  into  the  form  of  the  free- 
space  Green's  function  for  an  isospeed  medium  given  by 

-jkR 
<P(D  =  ^r.  (4.25) 

In  this  chapter,  we  will  refer  to  Eq.  (4.25)  as  the  exact  solution  for  the  free-space,  isospeed 
medium  since  no  simplifying  approximations  were  made  to  obtain  this  result.  Examination 
of  Eq.  (4.25)  shows  that  under  these  circumstances,  the  propagating  wave  undergoes 
spherical  spreading,  i.e.,  the  amplitude  of  the  acoustic  field  is  inversely  proportional  to 
radial  distance  between  the  transmitter  and  the  field  point.  This  radial  distance  is 
commonly  referred  to  as  the  Range-Line-Of-Sight  (RLOS). 

Reviewing  the  development  leading  up  to  the  expression  for  the  overall  system  complex 
frequency  response  due  to  a  unit  amplitude  point  source  given  by 

H(f,r)  =  27t    f     HM(f,fr,yo;y)J0(27tfrrm)e'j27tfY(yo)(y"yn)frdfr,  (4.33) 

remember  that  it  is  mathematically  the  same  as  the  acoustic  field  given  by  Eq.  (4.25).  As  a 
result,  the  overall  system  complex  frequency  response  should  behave  the  same  as  the 
acoustic  field  due  to  a  point  source.  Accordingly,  if  we  plot  RLOS  versus  the  magnitude 
of  the  overall  system  complex  frequency  response  on  log-log  axes  for  a  particular 
frequency  component  of  the  transmitted  signal,  we  would  expect  the  result  to  be  a  straight 


line.  Further,  due  to  the  isospeed  nature  of  the  medium,  we  would  expect  all  propagating 
frequency  components  to  behave  the  same  since  we  are  only  concerned  in  this  simulation 
with  the  effects  of  spreading  and  have  ignored  losses  due  to  absorption. 

B.   RESULTS  OBTAINED  BY  IGNORING  EVANESCENT  WAVES 

An  initial  investigation  was  carried  out  for  an  isospeed  medium  with  a  sound  speed  of 
1475  m/sec.  For  a  transmit  signal,  we  used  the  waveform  of  Fig.  5.1b,  which  was  a 
Lanczos  smoothed,  rectangular  amplitude  weighted  CW  pulse  with  the  following 
parameters : 

•  AMP  =  1.0,  TP  =  12.5  msec,  PRF  =  40.0  Hz,  FC  =  400.0  Hz 

•  minimum  frequency  =  160.0  Hz,  maximum  frequency  =  640  Hz. 

For  this  scenario,  the  longest  wavelength  (which  corresponds  to  the  lowest  frequency 
component  present)  is  equal  to  9.219  m,  while  the  shortest  wavelength  (which  corresponds 
to  the  highest  frequency  component  present)  is  equal  to  2.305  m.  Since  the  evanescent 
waves  are  exponentially  decaying  instead  of  being  oscillatory,  they  are  considered  as  being 
non-propagating  and  are  usually  ignored  at  distances  greater  than  a  couple  of  wavelengths 
from  their  source.  For  the  case  just  described,  we  would  expect  evanescent  waves  to  be 
insignificant  at  ranges  greater  than  say  fifty  meters.  Applying  this  assumption  to  our 
expression  for  the  overall  system  complex  frequency  response,  i.e.,  ignoring  the 
integration  over  the  region  of  evanescent  waves,  results  in 

/•f/co 
H(f,r)  =  27t    I  HM(f,fr,yo;y)J0(27tfrrm)e-j27rfY(yoXy-yn)frdfr,  (6.3) 
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which  is  an  approximation  of  the  full-wave  solution  since  the  region  of  integration  is  no 
longer  infinite.  We  will  refer  to  any  overall  system  complex  frequency  response  that  has 
such  a  limited  region  of  integration  as  our  approximate  (full-wave)  solution  to  the  ocean 
medium  propagation  problem. 

It  should  be  stressed  at  this  point  that  the  computer  implementation  of  the  overall 
system  complex  frequency  response  has  been  kept  as  general  as  possible.  This  means  that 
instead  of  having  to  program  separate  forms  of  Eq.  (3.14)  for  different  scenarios,  e.g.,  Eq. 
(6.3),  just  the  generalized  form  of  the  overall  system  complex  frequency  response  given  by 
Eq.  (3.14)  is  implemented  in  the  computer  code.  This  generalized  form  can  then  be 
modified  as  required,  by  a  number  of  user  defined  variables.  To  obtain  the  overall  system 
complex  frequency  response  given  by  Eq.  (6.3),  we  merely  set  the  user  defined  variables 

•  transmit  array  :  MT  =  1 ,  NT  =  1 ,  STEER  =  FALSE 

•  integration  routine  :  RATIO  =1.0 

•  sound  speed  profile  :  NGRAD  =  0 

where  MT  and  NT  describe  the  dimensions  (number  of  elements)  of  the  transmit  array,  as 
discussed  in  Chapter  2,  and  STEER  is  a  logical  variable  that  designates  whether  any 
beamsteering  is  to  be  done.  If  no  beamsteering  is  to  be  done,  then  the  transmit  array 
complex  weights,  given  by  Eq.  (2.25),  are  set  equal  to  unity.  The  variable,  RATIO, 
defined  as 

RATIO  =  frc(y0)/f  ,  (6.4) 
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determines  the  upper  limit  of  the  integration  with  respect  to  fr ,  while  the  variable  NGRAD 
set  equal  to  zero  designates  an  isospeed  medium  and  the  program  uses  the  closed  form 
expression  for  the  phase  integral  appearing  in  Eq.  (3.14). 

To  demonstrate  the  validity  of  the  preceding  discussion,  we  need  to  determine  if  our 
derived  overall  system  complex  frequency  response  does  in  fact  reduce  to  the  isospeed, 
free-space,  Green's  function  and  show  the  expected  (1  /  RLOS)  fall-off  in  the  magnitude  of 
the  acoustic  field.  In  other  words,  we  need  to  compare  the  results  given  by  the 
approximate  full-wave  and  exact  solutions.  Accordingly,  the  magnitude  of  the  overall 
system  complex  frequency  response  versus  range-line-of-sight  is  plotted  in  Fig.  6.1  for 
three  different  frequencies.  The  frequencies  chosen  correspond  to  the  minimum,  carrier, 
and  maximum  frequencies  of  the  transmitted  signal.  Additionally,  the  approximate  full- 
wave  and  exact  values  of  the  magnitude  of  the  acoustic  field  at  the  carrier  frequency  of  400 
Hz  and  for  the  values  of  RLOS  used  in  Fig.  6.1  are  detailed  later  in  Table  6.3. 

As  can  be  seen  from  Fig.  6.1,  the  magnitude  of  the  acoustic  field  has  been  determined 
for  values  of  RLOS  ranging  from  ten  meters  to  ten  thousand  meters.  This  is  equivalent  to 
a  range  of  one  to  thousands  of  wavelengths,  depending  on  which  of  the  frequency 
components  we  are  considering.  From  our  preceding  discussion,  we  would  expect  that 
the  acoustic  field  calculated  using  the  approximate  full-wave  solution,  given  by  Eq.  (6.3), 
will  not  reduce  exactly  to  the  isospeed,  free-space,  Green's  function  for  ranges  less  than 
about  fifty  meters.  This  is  because  we  have  assumed  that  evanescent  waves  have  a 
noticeable  effect  in  this  region,  which  Eq.  (6.3)  ignores.  Conversely,  would  we  expect 
close  agreement  for  ranges  greater  than  fifty  meters  as  we  have  assumed  that  the  effect  of 
evanescent  waves  in  this  region  is  negligible.  In  terms  of  the  log-log  plot  of  magnitude 
versus  RLOS  shown  in  Fig.  6.1,  we  would  expect  the  graph  to  be  a  straight  line  for  ranges 
greater  than  fifty  meters,  and  to  deviate  away  from  this  ideal  for  shorter  ranges. 
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RLOS  (m) 

Fig.  6.1   Magnitude  of  the  acoustic  field  versus  RLOS  when 
the  effects  of  evanescent  waves  are  ignored. 
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It  can  readily  be  observed  from  Fig.  6.1  that  for  long  ranges,  the  graph  does  indeed 
approximate  a  straight  line  with  the  expected  slope,  and  for  short  ranges,  it  deviates  from 
this  straight  line.  These  results  appear,  at  first  glance,  to  be  what  we  expected.  However, 
note  the  range  at  which  the  plotted  curve  starts  to  approximate  a  straight  line.  This  range  is 
about  five  hundred  meters,  or  ten  times  the  range  that  we  thought  would  be  significant. 
Alternatively,  this  range  represents  between  54  and  216  wavelengths,  depending  on  which 
of  the  three  frequency  components  we  are  considering.  This  would  suggest  that  the  effects 
of  evanescent  waves  are  indeed  still  significant  at  ranges  greater  than  a  couple  of 
wavelengths  from  their  source. 

Consider  also,  the  actual  values  of  the  approximate  acoustic  field  that  are  contained  in 
Table  6.3.  Even  at  a  range  of  ten  thousand  meters,  the  approximate  full-wave  and  exact 
solutions  calculated  for  the  carrier  frequency,  still  differ  by  almost  1%.  For  the  carrier 
frequency,  this  range  corresponds  to  2,712  wavelengths,  which  is  well  beyond  where  we 
expected  evanescent  waves  to  have  an  effect. 

Close  inspection  of  Fig.  6.1  does  indicate,  however,  that  the  low  frequency  term,  i.e., 
160  Hz,  is  the  one  that  is  behaving  the  most  erratically.  This  is  what  we  would  expect 
with  regards  to  the  deviation  from  the  straight  line  predicted  by  the  exact  solution  since  this 
frequency  corresponds  to  the  longest  wavelength.  Consequently,  for  any  given  distance, 
this  is  the  evanescent  component  that  would  have  decayed  the  least.  Since  we  are  ignoring 
evanescent  waves  in  this  case,  this  is  the  frequency  component  that  should  deviate  most 
from  the  exact  result.  Based  on  these  observations  and  trends,  we  will  assume  that  our 
derivation  of  the  overall  system  complex  frequency  response  has  not  been  proved  in  error 
so  far,  and  that  it  is  the  assumption  regarding  the  significance  of  evanescent  waves  that  is 
incorrect. 
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Even  though  our  approximate  solution  results  deviate  from  those  predicted  by  the  exact 
solution,  is  the  amount  of  deviation  really  significant  ?  Remember  that  one  of  the  outputs 
that  we  wanted  from  the  computer  simulation  was  the  time-domain  electrical  signal  at  each 
element  in  the  receive  array,  which  would  be  used  to  show  pulse  distortion  due  to  the 
effects  of  dispersion,  and  as  input  data  to  space-time  signal  processing  algorithms.  To 
demonstrate  this  pulse  distortion,  we  need  to  compare  the  original  transmitted  electrical 
signal  with  the  received  electrical  signal  and  observe  the  changes  that  propagation  through 
the  ocean  medium  has  caused.  Such  a  comparison  can  be  made  using  Figs.  6.2a  and 
6.2b.  Note  that  6.2b  is  the  received  electrical  signal  at  a  RLOS  of  100  m.  Referring  to 
Fig.  6.1,  we  see  that  this  value  of  RLOS  is  one  for  which  the  calculated  acoustic  field 
varies  significantly  with  frequency.  This  particular  RLOS  was  chosen  to  represent  a  near 
worst  case  situation  and  will  demonstrate  the  significance,  if  any,  of  the  inexactness  of  the 
approximate  full-wave  solution  when  evanescent  waves  are  ignored. 

It  can  be  seen  from  comparing  Figs.  6.2a  and  6.2b,  that  the  output  electrical  signal  is 
indeed  significantly  different  from  the  transmitted  signal.  Recall  that  there  should  be  no 
dispersion  in  an  isospeed,  free-space  medium.  This  means  that  the  received  electrical 
signal  should  be  a  scaled  replica  of  the  transmitted  electrical  signal,  i.e.,  identical  in  pulse 
shape.  The  problem  with  the  approximate  full-wave  solution  that  produced  the  output 
electrical  signal  shown  in  Fig.  6.2b  is  that,  if  we  repeated  the  problem  for  a  non-isospeed 
medium  in  which  we  do  expect  pulse  distortion  due  to  dispersion,  we  would  not  know 
how  much  of  the  distortion  was  due  to  the  dispersion  and  how  much  was  due  to  the 
approximate  nature  of  our  solution. 

Next,  consider  the  previous  problem  repeated  for  a  RLOS  of  1,000  m.  From  Fig.  6.1 
it  can  be  seen  that  for  this  range,  the  acoustic  field  behaves  a  lot  more  as  predicted  and  is 
less  dependent  on  frequency.    This  means  that  there  is  much  closer  agreement  between  the 
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Fig.  6.2a      Transmitted  electrical  signal. 
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Fig.  6.2b      Received  electrical  signal,  RLOS  =  100  in,  when  the 
effects  of  evanescent  waves  are  ignored. 


approximate  full-wave  and  exact  solutions  than  existed  for  RLOS  =  100  m.  The  output 
electrical  signal  for  this  case  is  shown  in  Fig  6.2c,  where  the  transmitted  electrical  signal  is 
the  same  as  for  the  previous  case  shown  in  Fig.  6.2b.  As  expected,  the  output  electrical 
signal  indeed  looks  a  lot  more  like  the  transmitted  electrical  signal  than  was  the  case  in  Fig. 
6.2a.  In  fact  the  amount  of  distortion  can  be  considered  imperceptible  for  RLOS  =  1,000 
m,  and,  for  all  practical  purposes,  Fig.  6.2c  is  a  scaled  version  of  die  transmitted  signal. 

Based  on  these  results  for  the  received  time-domain  electrical  signal,  it  appears  that  the 
approximate  full-wave  solution  is  adequate  for  demonstrating  the  effects  of  pulse  distortion 
for  signals  containing  no  frequency  components  below  160  Hz,  where  the  receive  element 
is  at  a  range  greater  than  1  ,(KX)  m  from  the  source.  This  situation  should  be  acceptable  for 
long  range  propagation  problems,  but  for  any  short  range  calculations,  it  is  obviously 
inadequate. 
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Fig.  6.2c     Received  electrical  signal,  RLOS  =  1,000  rn,  when 
the  effects  of  evanescent  waves  are  ignored. 
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Consider  next,  the  other  use  of  the  time-domain  electrical  signal  at  each  element  in  the 
receive  array,  that  is,  as  input  data  to  space-time  signal  processing  algorithms.  As  stated 
earlier,  the  computer  code  developed  for  this  ocean  medium  propagation  problem  does  not 
perform  any  signal  processing  on  the  received  electrical  signal  at  each  point  in  the  array,  but 
rather,  the  data  is  written  to  a  file  so  that  it  can  be  operated  on  by  ancillary  programs.  As  a 
further  validation  of  the  generated  results,  we  can  operate  on  the  time-domain  data  with  a 
beamforming  algorithm  that  localizes  the  position  of  the  source,  and  compare  the  computed 
and  actual  source  position.  The  algorithm  used  was  an  adaptive  beamforming  and  non- 
linear estimation  algorithm  developed  by  Breckenridge  [Ref.  16]  that  localizes  a  source  in 
range,  depression  angle,  and  azimuthal  angle,  according  to  the  geometry  outlined  in  Fig. 
6.3.  From  Fig.  6.3  it  should  be  noted  that  the  localization  of  the  source,  calculated  by  this 
adaptive  algorithm,  is  relative  to  the  center  element  in  the  receive  array. 

The  data  generated  from  the  two  cases,  RLOS  =  100  m  and  RLOS  =  1,000  m,  was 
processed  with  this  algorithm  using  the  output  electrical  signal  from  each  of  the  elements  in 
a  3x3  receive  array,  and  the  results  are  shown  in  Table  6.1.  The  values  RLOS  (X)  and 
RLOS  (Y)  detailed  in  the  table  correspond  to  estimates  of  the  RLOS  calculated  from 
information  along  the  X  and  Y  axes,  respectively.  From  the  results  contained  in  Table 
6.1,  we  can  see  that  the  estimated  values  from  the  beamforming  algorithm  are  closer  to  the 
true  values  for  the  case  RLOS  =  1,000  m  than  for  the  case  RLOS  =  100  m,  which  is 
consistent  with  our  previous  observations  in  this  section. 

An  interesting  feature  of  the  results  is  the  way  that  RLOS  00  differs  radically  from  the 
true  value  of  RLOS,  while  RLOS  (X)  is  fairly  close  to  the  true  value.  Note  that  this  occurs 
for  both  values  of  RLOS  considered,  that  is,  RLOS  (Y)  does  not  improve  as  the  value  of 
RLOS  is  increased.  This  phenomena  is  in  contrast  to  the  trends  observed  earlier  in  this 
section.    The  fact  that  the  effects  of  evanescent  waves  are  being  ignored  means  that  we  are, 
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Fig.  6.3  Geometry  for  adaptive  beamforming  algorithm. 


RLOS 
(m) 

true 

RLOS  (X) 
(m) 
est 

RLOS  (Y) 
(m) 

est 

0 

(deg) 

true 

0 

(deg) 
est 

(deg) 
true 

4> 

(deg) 
est 

100 
1000 

115.02 
960.71 

29.54 
284.619 

30.000 
30.000 

28.274 
29.942 

270.00 
270.00 

270.00 
270.00 

Table  6.1   Adaptive  beamforming  results  when 
evanescent  waves  are  ignored. 


in  essence  affecting  the  wave  solution  in  the  Y  direction  to  a  much  greater  degree  than  in  the 
X  direction.  Recall  that  our  form  of  the  overall  system  complex  frequency  response,  given 
by  Eq.  (3.14),  was  developed  assuming  a  sound-speed  profile  that  was  a  function  of 
depth.  In  our  development  we  showed  that  this  meant  Eq.  (3.14)  was  axisymmetric.  If 
the  overall  system  complex  frequency  response  is  axisymmetric,  then  the  relative  behavior 
of  the  acoustic  field  in  the  X  direction  should  be  minimally  affected  by  the  behavior  of  the 
acoustic  field  in  the  Y  direction.  Consequently,  if  there  is  a  problem  with  the  description 
of  the  full-wave  solution  in  the  Y  direction,  as  in  this  case,  the  value  of  RLOS  (X)  should 
still  be  fairly  accurate,  which  is  the  result  that  we  observed. 

These  results  also  demonstrate  that  while  the  time-domain  output  electrical  signal  from  a 
single  receive  array  element  may  look  correct  to  the  eye  and  suggest  that  the  accuracy  of  the 
generated  data  is  adequate,  it  may  in  fact  not  be  accurate  or  correct  enough  for  signal 
processing  purposes.  This  suggests  that  observing  the  output  electrical  signal  is  not  a 
good  way  of  validating  the  performance  of  the  approximate  full-wave  solution. 

The  remaining  output  that  we  wanted  from  this  computer  simulation  was  the  magnitude 
and  phase  of  the  complex  acoustic  field  incident  upon  a  receive  array  as  a  function  of 
frequency  and  spatial  coordinates.  A  convenient  way  of  viewing  these  results,  apart  from 
a  three-dimensional  plot,  is  to  tabulate  the  magnitude  and  phase  of  the  acoustic  field  in  the 
same  order  as  the  elements  of  the  receive  array.  Consider  a  3x3  element  receive  array. 
We  will  adopt  the  element  designation  scheme  shown  in  Fig.  6.4  to  number  the  elements 
(  m,n  )  contained  in  the  array. 

The  values  of  the  magnitude  and  phase  (in  degrees)  of  the  acoustic  field  are  shown 
in  Tables  6.2a  and  6.2b  for  the  cases  RLOS  =  100  m  and  RLOS  =  1,000  m,  respectively. 
Note  that  the  first  value  is  the  magnitude  and  the  second  value  is  the  phase.  Given  the 
placement  of  the  receive  array  relative  to  the  source,  i.e.,  XR  =  0  and  YR  =  1025  m,  which 
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Fig.  6.4   Element  designation  for  a  3x3  receive  array. 

is  below  the  source,  we  would  expect  the  magnitude  of  the  acoustic  field  to  decrease  with 
increasing  n.  This  corresponds  to  reading  down  the  table,  or  equivalently,  increasing  the 
depth  of  the  receive  element  or  field  point.  We  would  also  expect  the  value  of  the  acoustic 
field  to  decrease  with  increasing  absolute  value  of  m.  This  corresponds  to  reading  across 
the  table  out  from  the  center,  or  equivalently,  increasing  the  displacement  in  the  X  direction 
of  the  element  or  field  point.  We  expect  the  magnitude  of  the  acoustic  field  to  decrease  in 
both  of  these  situations  since  their  effect  is  to  increase  the  radial  range  between  the  source 
and  the  receive  element ,  or  field  point. 

From  these  two  tables  it  can  be  seen  that  all  the  values  for  the  acoustic  field  do  not 
behave  as  expected.  In  Table  6.2a,  the  first  two  rows  of  results  have  the  magnitude  of  the 
acoustic  field  increasing  as  we  move  away  from  m  =  0.  In  Table  6.2b  on  the  other  hand,  it 
is  the  third  or  final  row  of  values  that  behaves  in  this  fashion.  In  both  cases  however,  the 
magnitude  of  the  acoustic  field  decreases  as  the  depth  is  increased,  as  expected.    As  for  the 
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n/m 

-1 
0 

1 


8.69822  e-4 
-168.608 


8.47479  e-4 
139.860 


7.82107  e-4 
84.800 


8.68689  e-4 
-167.994 


8.47246  e-4 
140.423 


7.82997  e-4 
85.376 


8.69822  e-4 
-168.608 


8.47479  e-4 
139.860 


7.82107  e-4 
84.800 


Table  6.2a     Acoustic  field  at  array  elements,  RLOS  =  100  m. 


n/m 

1 

0 

-1 

-1 

8.00331  e-5 
169.065 

8.00375  e-5 
169.130 

8.00331  e-5 
169.065 

0 

7.97398  e-5 
113.082 

7.97420  e-5 
113.150 

7.97398  e-5 
113.082 

1 

7.93029  e-5 
56.758 

7.93010  e-5 
56.826 

7.93029  e-5 

56.758 

Table  6.2b    Acoustic  field  at  array  elements,  RLOS  =  1,000  m. 
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results  obtained  from  the  use  of  the  adaptive  beamforming  algorithm,  observing  the  values 
of  the  complex  acoustic  field  as  a  function  of  frequency  and  spatial  coordinates  alerts  us  to 
the  fact  that  the  data  produced  by  our  approximate  full-wave  solution  is  inaccurate. 

From  the  results  of  this  section,  we  can  formulate  two  primary  conclusions.  First,  the 
assumption  that  the  effects  of  evanescent  waves  past  a  distance  of  a  couple  wavelengths  is 
insignificant  is  an  invalid  one,  at  least  for  the  situation  we  are  modelling  and  for  the  full- 
wave  method  used.  Second,  comparing  the  transmitted  and  received  electrical  signals  for 
the  free-space,  isospeed  medium  case  is  not  sufficient  to  determine  whether  or  not  the 
approximate  full-wave  solution,  is  accurate.  To  gauge  the  performance  of  the  approximate 
solution  we  need  to  either  observe  the  acoustic  field  as  a  function  of  spatial  coordinates  for 
an  appreciation  for  whether  the  solution  is  behaving  correctly,  or  apply  signal  processing 
techniques  to  obtain  some  quantitative  measure  for  the  accuracy  of  the  solution. 

C.    RESULTS  OBTAINED  BY  INCLUDING  EVANESCENT  WAVES 

To  include  the  effects  of  evanescent  waves,  we  need  to  extend  the  region  of  integration 
with  respect  to  fr  in  Eq.  (6.3).  As  was  previously  discussed,  the  region  of  evanescent 
waves  extends  from  fr  equal  to  (f  /  c0)  to  infinity.  But  is  it  necessary  to  integrate  over  this 
entire  region  ?  After  some  trial  and  error,  it  was  found  that  integrating  with  respect  to  fr  up 
to  a  limit  of  1 .2  (f  /  c0)  was  adequate  to  give  very  good  results  down  to  a  RLOS  of  10  m. 
Recall  that  evanescent  waves  have  more  effect  at  short  ranges.  Therefore,  if  this  upper  limit 
of  integration  is  sufficient  for  a  RLOS  of  10  m,  then  it  will  also  be  sufficient  for  longer 
ranges.  Since  no  practical  reason  could  be  thought  of  as  to  why  we  would  be  interested  in 
ranges  less  than  10  m,  this  is  the  upper  limit  of  integration  that  we  adopted  and  is  that  used 
to  obtain  the  results  of  this  section.  Therefore,  our  approximate  full- wave  solution  is  now 
given  by 
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-1.2(f/c0) 
H(f,r)  =  2n   I  HM  (f,fr,y0;y)  J0(2idTTm )  e^^0^  '  Yn)fr  dfr.    (6.5) 


Recall  from  our  previous  discussion  that  the  region  of  integration  in  Eq.  (6.5)  is,  for 
computational  purposes,  divided  into  two  regions  to  account  for  the  change  in  the  form  of 
the  propagation  vector  in  the  Y  direction,  k  Y  (y),  when  fr  changes  from  being  less  than  to 
greater  than  f  /  c0 .  This  empirically  determined  value  for  the  upper  limit  of  integration 
agrees  with  observations  made  by  Schmidt  [Ref.  7:p.  31]  with  regard  to  wavenumber 
integrations. 

For  the  analysis  contained  in  this  section,  the  scenario  is  the  same  as  for  the  previous 
case  where  we  ignored  the  effects  of  evanescent  waves,  i.e.,  same  transmit  signal,  array 
geometry  and  ocean  medium.  To  allow  for  this  extended  region  of  integration,  all  that  has 
to  be  done  is  to  set  the  variable  RATIO  equal  to  1 .2  instead  of  1 .0.  This  demonstrates  one 
of  the  advantages  in  having  kept  our  initial  analysis  and  computer  code  implementation  as 
general  as  possible.  It  is  a  very  simple  matter  to  alter  the  problem  that  requires  no 
programming  changes. 

The  magnitude  of  the  overall  system  complex  frequency  response  was  determined  for  a 
range  of  values  for  RLOS  and  for  three  frequencies  corresponding  to  the  minimum,  carrier, 
and  maximum  frequency  components  of  the  transmitted  signal.  It  was  found,  however, 
that,  unlike  the  previous  case,  there  was  no  discemable  frequency  dependence  in  the  values 
for  the  magnitude  of  the  overall  complex  frequency  response  determined  from  the 
approximate  full-wave  solution  given  by  Eq.  (6.5).  As  was  previously  discussed,  this 
frequency  independence  is  the  result  predicted  from  the  exact  form  of  the  isospeed,  free- 
space  Green's  function.  Since  the  results  for  the  three  frequency  components  were 
identical,  only  the  results  for  the  carrier  frequency,  400  Hz,  are  shown  in  the  log-log  plot 


96 


of  acoustic  field  versus  RLOS  shown  in  Fig.  6.5.  Additionally,  the  approximate  full-wave 
and  exact  values  of  the  magnitude  of  the  acoustic  field  for  this  frequency  and  values  of 
RLOS  used  in  Fig.  6.5  are  detailed  in  Table  6.3. 

From  Fig.  6.5,  it  can  be  observed  that  the  log-log  plot  of  the  acoustic  field  magnitude 
versus  RLOS  is  identical  in  slope  to  the  straight  line  predicted  by  the  isospeed,  free-space, 
Green's  function,  over  the  entire  range  of  values  of  RLOS  considered.  Comparing  Figs. 
6.1  and  6.5,  it  can  easily  be  observed  that  there  has  been  a  radical  improvement  due  to  the 
inclusion  of  a  small  region  of  evanescent  waves.  The  results  from  Fig.  6.5  suggest  that 
the  upper  limit  of  integration  chosen  in  the  approximate  full-wave  solution  given  by  Eq. 
(6.5)  is  adequate  for  values  of  RLOS  down  to  at  least  10  m. 

Inspection  of  Table  6.3  shows  that  with  the  inclusion  of  this  region  of  evanescent 
waves,  the  magnitude  of  the  acoustic  field  given  by  the  approximate  full-wave  solution  in 
fact  agrees  with  the  exact  values  predicted  by  the  isospeed,  free-space,  Green's  function  to 
within  five  significant  figures.  The  comparison  was  limited  to  this  number  of  significant 
figures  since  the  numerical  integration  was  only  carried  out  to  sufficient  terms  to  give  a 
maximum  of  a  0.01%  relative  error  in  the  value  of  the  integral.  Alternatively  stated,  the 
value  of  the  acoustic  field  calculated  from  Eq.  (6.5)  is  accurate  to  within  ±0.01%. 

Let  us  now  consider  the  displaying  of  the  received  electrical  signal  for  the  purposes  of 
showing  pulse  distortion  due  to  dispersion.  Comparing  Figs.  6.6a  and  6.6b,  it  can  be 
observed  that  there  is  no  discernable  pulse  distortion,  which  is  what  theory  predicts  for 
wave  propagation  in  a  free-space,  isospeed  medium.  This  result  was  also  expected  from 
experience  with  the  analysis  of  the  previous  case,  in  which  we  ignored  the  effects  of 
evanescent  waves.  In  that  case,  even  though  there  was  some  frequency  dependence  in  the 
value  of  the  acoustic  field  calculated  using  the  approximate  full-wave  solution  of  Eq.  (6.3), 
there  was  very  little  discernable  distortion  in  the  received  electrical  signal.    For  the  acoustic 
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RLOS  (m) 

Fig.  6.5  Magnitude  of  the  Acoustic  field  versus  RLOS  when 
the  effects  of  evanescent  waves  are  included. 
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RLOS(m) 

1/4tiR 

1  H(f,r)  1 

1  H(f,r)  1 

( no  evanescent 
waves) 

( incl.  evanescent 
waves ) 

10 

7.9577  e-3 

7.9702  e-3 

7.9577  e-3 

50 

1.5915  e-3 

1.3983  e-3 

1.5915  e-3 

100 

7.9577  e-4 

8.4725  e-4 

7.9577  e-4 

200 

3.9789  e-4 

3.9676  e-4 

3.9789  e-4 

300 

2.6526  e-4 

2.5923  e-4 

2.6526  e-4 

500 

1.5915  e-4 

1.5663  e-4 

1.5915  e-4 

1,000 

7.9577  e-5 

7.9742  e-5 

7.9577  e-5 

2,000 

3.9789  e-5 

3.9091  e-5 

3.9789  e-5 

3,000 

2.6526  e-5 

2.5978  e-5 

2.6526  e-5 

5,000. 

1.5915  e-5 

1.5753  e-5 

1.5915  e-5 

10,000 

7.9577  e-6 

8.0182  e-6 

7.9577  e-6 

Table  6.3   Comparison  between  exact  and  predicted 
magnitude  values  of  the  acoustic  field. 
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field  calculated  using  the  more  accurate  approximate  full-wave  solution  given  by  Eq.  (6.5), 
we  would  expect  our  received  electrical  signal  to  look  even  more  like  a  scaled  version  of  the 
transmit  signal.  This  was,  in  fact,  the  observed  result. 

The  beamforming  algorithm  was  then  applied  to  the  output  electrical  signals  from  a  3x3 
receive  array,  as  in  the  previous  case,  with  the  results  being  summarized  in  Table  6.4.  The 
first  thing  to  notice  about  these  results  is  that  unlike  the  previous  case  that  ignored 
evanescent  waves,  both  RLOS  (X)  and  RLOS  (Y)  give  accurate  values  for  the  radial 
distance  between  the  source  and  the  center  of  the  receive  array.  Additionally,  the  estimated 
position  angles  are  very  accurate. 

Finally,  an  inspection  of  the  complex  acoustic  field  at  the  elements  of  the  receive  array 
is  shown  in  Table  6.5.  It  can  be  clearly  seen  that  the  values  of  the  magnitude  conform  to 
the  expected  pattern,  as  described  for  the  previous  section.  Note  that  the  change  in 
magnitude,  reading  down  the  table,  is  much  greater  than  when  reading  across  the  table. 
While  this  may  appear  to  be  an  error  at  first,  it  is  in  fact  correct,  and  is  due  to  the 
positioning  of  the  array.  Recall  that  the  receive  array  is  positioned  such  that  XR  =  0, 
YR  =  1050  m,  i.e.,  the  receive  array  is  50  m  below  the  transmit  array,  and  RLOS  =  100  m. 
Also  recall  that  the  distance  between  elements  in  the  X  and  Y  directions  is  the  same. 
Accordingly,  it  can  be  easily  shown  that  the  change  in  radial  distance  between  the  source 
and  a  receive  element  is  much  larger  for  an  incremental  change  in  the  Y  direction  than  for 
the  same  change  in  the  X  direction.  As  a  result,  the  change  in  the  magnitude  of  the 
acoustic  field  will  be  similarly  affected. 

From  the  results  of  this  chapter  so  far,  we  have  seen  that  it  is  not  a  valid  assumption  to 
ignore  the  effects  of  evanescent  waves  when  developing  an  approximation  of  the  full- wave 
solution  to  an  ocean  medium  propagation  problem.  However,  it  appears  that  we  only  have 
to  include  a  small  portion  of  the  evanescent  region  for  an  approximate  full-wave  solution  to 
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Fig.  6.6a      Transmitted  electrical  signal. 
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RLOS:  100.0  M   BLOS:  80. 0  DEC 


Fig.  6.6b      Received  electrical  signal,  RLOS  =  100  m,  when  the 
effects  of  evanescent  waves  are  included/ 
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be  valid.  From  these  results  we  will  further  assume  that  the  region  of  integration  with 
respect  to  fr  ,  given  by 

0<fr<1.2(f/c0),  (6.6) 

for  the  overall  system  complex  frequency  response,  given  by  Eq.  (3.14),  is  adequate  for 
sound- speed  profiles  other  than  the  isospeed  case  investigated  in  this  chapter.  We  will 
utilize  this  assumption  in  the  following  chapter,  when  we  investigate  the  behavior  of  the 
acoustic  field  in  a  free-space  medium  characterized  by  a  linear  sound- speed  profile  with  a 
single  gradient.  We  have  also  seen  that  the  accuracy  of  the  approximate  full-wave  solution 
is  a  function  of  the  numerical  integration  used  to  calculate  the  overall  system  complex 
frequency  response.  There  is  a  trade-off  between  computation  time  and  accuracy  that  will 
be  addressed  in  the  final  section  of  this  chapter. 

D.   PLOTTING  OF  THE  3-D  COMPLEX  ACOUSTIC  FIELD 

So  far  in  this  chapter,  all  our  validation  has  taken  place  using  a  similar  geometry,  i.e., 
XR  =  0  and  YR  >  YT.  This  means  that  the  receiver  has  always  been  deeper  than  the 
transmitter,  and  has  had  no  offset  in  the  X  direction.  To  have  any  confidence  that  the 
theory  and  computer  implementation  of  the  overall  system  complex  frequency  response  is 
correct,  we  need  results  for  every  permutation  of  the  X  and  Y  direction  offsets  of  the 
receiver  with  respect  to  the  source.  In  other  words,  we  need  to  place  the  receiver  in  a 
variety  of  positions,  e.g.,  positive  X  offset  and  above  the  source,  to  fully  test  the 
algorithm. 

Since  one  of  the  required  outputs  is  the  magnitude  and  phase  of  the  complex  acoustic 
field  at  the  locations  of  the  array  elements  for  use  by  matched-field  processing  algorithms,  a 
3-D  plotting  routine  was  implemented  and  added  as  a  visual  aid  to  help  conceptualize  what 
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RLOS 
(m) 
true 

RLOS(X) 
(m) 
est 

RLOS(Y) 
(m) 

est 

e 

(deg) 
true 

e 

(deg) 
est 

(deg) 
true 

(deg) 
est 

100 

99.984 

99.977 

30.000 

30.004 

270.00 

270.00 

Table  6.4   Adaptive  beamforming  results  when 
evanescent  waves  are  included. 


n/m 

1 

0 

-1 

-1 

8.00292  e-4 
-167.604 

8.00346  e-4 
-166.952 

8.00191  e-4 
-167.604 

0 

7.95722  e-4 
136.639 

7.95774  e-4 
137.287 

7.95722  e-4 
136.639 

1 

7.91125  e-4 
79.909 

7.91 177  e-4 
80.554 

7.91125  e-4 
79.909 

Table  6.5    Acoustic  field  at  array  elements,  RLOS  =  100  m. 
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the  acoustic  field  looked  like.  This  routine  plots,  on  separate  graphs,  the  magnitude  and 
phase  of  the  overall  system  complex  frequency  response  for  a  particular  frequency  versus 
the  relative  position  of  each  of  the  array  elements.  Recall  from  our  earlier  development  that 
the  value  of  the  frequency  spectrum  of  the  acoustic  field  at  any  point  is  just  the  value  of  the 
overall  system  complex  frequency  response  times  the  frequency  spectrum  of  the  transmitted 
electrical  signal.  By  then  looking  at  the  plots,  particularly  those  of  the  magnitude  of  the 
acoustic  field,  it  can  easily  be  determined  if  the  field  exhibits  the  appropriate  orientation  and 
curvature  for  the  given  relative  positions  of  the  source  and  the  receive  array. 

Note  that  the  planar  array  of  point  sources  we  are  calling  the  transmitter  does  not 
necessarily  have  to  be  a  physical  transducer  array.  It  could  represent  sampled  values  of  the 
acoustic  field  at  some  spatial  location  (or  grid  of  locations)  due  to  some  other  source.  In 
this  situation,  the  array  could  be  thought  of  as  detecting  an  acoustic  signal  and  reradiating 
that  acoustic  signal  without  alteration. 

As  was  previously  discussed,  we  need  to  generate  a  series  of  plots  of  the  complex 
acoustic  field  for  various  combinations  of  source  /  receiver  locations.  A  series  of  plots 
were  generated  for  the  six  different  permutations  of  the  X  and  Y  offsets  of  the  receive  array 
relative  to  the  source.  For  all  these  plots,  the  source  was  located  at  a  depth  of  1,000  m. 
Recall  that  the  transmit  array  (or  single  point  source  in  this  case)  is  centered  within  the 
coordinate  axis  system  such  that  XT  and  ZT  are  both  equal  to  zero  and  YT  is  the  depth 
below  the  sea  surface  (see  Fig.  2.3).  The  receiver  was  a  9  x  9  planar  array  of  point 
sources,  which  resulted  in  the  computation  of  the  overall  system  complex  frequency 
response  for  81  different  spatial  locations.  These  values  of  the  overall  system  complex 
frequency  response  were  calculated  for  a  frequency  of  400  Hz,  which  is  the  carrier 
frequency  of  the  transmitted  electrical  signal  shown  in  Figs.  6.2a  and  6.6a.  The  RLOS, 
i.e.,  the  radial  distance  between  the  source  point  and  the  center  of  the  receive  array,  was 
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assumed  to  be  only  50  m  in  order  to  reduce  computer  processing  time.  This  short  range 
was  chosen  so  that  we  could  process  data  a  reasonable  sized  array  which  would  allow  us 
to  observe  curvature  in  the  acoustic  field,  i.e.,  so  that  the  plots  would  not  just  look  like  flat 
planes.  The  array  size  was  chosen  as  being  adequate  to  produce  a  3-D  plot,  for  the  RLOS 
as  previously  chosen,  in  which  the  orientation  and  curvature  of  the  acoustic  field  values 
was  easily  discernable.  A  constant  speed  of  sound  of  1475  m  /  sec  was  assumed  for  all 
cases.  In  the  following  3-D  plots,  the  X  and  Y  offsets  AY  =  YR  -  YT  and  AX  =  XR  are 
used  to  describe  the  relative  positions  of  the  centers  of  the  transmit  and  receive  arrays. 

A  series  of  six  sets  of  plots  were  generated,  showing  both  the  magnitude  and  phase  of 
the  overall  system  complex  frequency  response,  for  the  following  positions  of  the  receive 
array  : 

•  Fig.  6.7  YR  =  1000  m,    XR  =  0  m,  (AY  =  0,  AX  =  0) 

•  Fig.  6.8  YR  =  1050  m,    XR  =  0  m,  (AY  =  50,  AX  =  0) 

•  Fig.  6.9  YR  =  1010  m,    XR  =  20  m,  (AY  =  10,  AX  =  20) 

•  Fig.  6.10  YR  =  980  m,    XR  =  10  m,  (AY  =  -  20,  AX  =  10) 

•  Fig.  6.1 1   YR  =  1020  m,  XR  =  -  20  m,  (AY  =  20,  AX  =  -  20) 

•  Fig.  6.12  YR  =  990  m,    XR  =  -  20  m,  (AY  =  -  10,  AX  =  -  20) 

Recall  that  the  parameters,  (XR,  YR,  ZR),  describe  the  position  of  the  center  element  of  the 
receive  array. 

Inspection  of  the  magnitude  plots  of  these  figures  clearly  show  that  the  computer  code 
does  indeed  handle  the  different  permutations  of  the  X  and  Y  offsets.  Consider  the 
magnitude  plot  shown  in  Fig.  6.10a.  Since  AY  is  negative  and  AX  is  positive,  then  from 
the  perspective  of  the  receiver  looking  towards  the  source,  the  receiver  is  located  above  and 
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CASE :  HOf-lOG   CW  FREQUENCY:    4Q0.0    HZ 

YT:     1OO0.0    M      XR:    0.0    M      YR:     1000.0    M      ZR:    50.0    M 

DXR^     1.152    M      DYR:     1.152    M 


Fig.  6.7a  Magnitude  of  the  overall  system  complex  frequency  response. 
(AY  =  0,AX  =  0) 
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CR5E:  MOMOG  Cl-J    FREQUENCY:  400.0  HZ 

YT:  1000.0  M   XR:  0.0  M   YR:  1000.0  M   ZR:  50.0  M 

DXRi-  1.152    11      OYR:     1.152    11 


Fig.  6.7b  Phase  of  the  overall  system  complex  frequency  response. 
(AY  =  0,AX  =  0) 


107 


CR5E:H0M0G   CW  FREQUENCY:    400.0    HZ 

YT:     1000.0    M      XR:    0.0    M      YR:     1050.0    M      ZR:    86.6    M 

DXR:'1.152    M      DYR:     1.152    M 


Fig.  6.8a  Magnitude  of  the  overall  system  complex  frequency  response. 
(AY  =  50,  AX  =  0) 
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CRSE: HOMOG   CW  FREQUENCY:    400.0    HZ 

YT:     1000.0    M      XR:    0.0    M      YR:     1050.0    M      ZR:    86.6    M 

DXR:     1.152    M      DYR:     1 . 152    M 


Fig.  6.8b  Phase  of  the  overall  system  complex  frequency  response. 
(AY  =  50,  AX  =  0) 
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CRSErHOMOG   CM  FREQUENCY:    400.0    HZ 

YT:     1000.0    M      XR:    20.0    M      YR:     1010.0    M      ZR:    44.7    M 

DXR:     1.152    M      DYR:     1.152    M 


Fig.  6.9a  Magnitude  of  the  overall  system  complex  frequency  response. 
(AY=  10,  AX  =  20) 
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CRSE:HOMOG  CW  FREQUENCY:    400.0    HZ 

YT:    1000.0    M      XR:    20.0    M      YR:     1010.0    M      ZR:    4A.7    M 

DXR:'   1.152    M      DYR:     1.  152    M 


Fig.  6.9b  Phase  of  the  overall  system  complex  frequency  response. 
(AY=  10,  AX  =20) 
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CnSE:HOMOG  CW    FREQUENCY:  400.0  HZ 

YT:  1000.0  M   XR:  10.0  M   YR:  980.0  M   ZR:  44.7  M 

DXRi- 1.152    M      OYR:     1.152    N 


Fig.  6.10a  Magnitude  of  the  overall  system  complex  frequency  response. 
(AY  =  -20,  AX  =  10) 


112 


CR5E : HOMOG   CW  FREQUENCY:    400.0    HZ 

YT:     1000.0    M      XR:     10.0    M      YR:    900.0    M      ZR:    44.7    M 

DXR:     1.152    M      DYR:     1.152    M 


Fig.  6.10b  Phase  of  the  overall  system  complex  frequency  response. 
(AY  =  -  20,  AX  =  10) 
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CnSErHOMOG   CW  FREQUENCY:    400.0    HZ 

YT:     1000.0    M      XR:-20.0    M      YR:     1020.0    M      ZR:    41.2    M 

DXR:'  1.152    M      DYR:     1.152    M 


Fig.  6.1  la  Magnitude  of  the  overall  system  complex  frequency  response. 
(AY  =  20,  AX  =  -  20) 
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Cn5E:H0M0G   CH  FREOUENCY:    400.0    HZ 

YT:    1000.0    M      XR:-20.0    M      YR:     1020.0    M      ZR:    41.2    M 

DXR:'  1.152    f1      DYR:     1  .152    M 


Fig.  6.1  lb  Phase  of  the  overall  system  complex  frequency  response. 
(AY  =  20,  AX  =  -20) 
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CnSE:HOMOG   CW  FREQUENCY:    400.0    HZ 

YT:     1000.0    II      XR:-20.0    M      YR:    990.0    M      ZR:    44.7    M 

DXR:'   1.152    M      DYR:     1.152    M 


Fig.  6.12a  Magnitude  of  the  overall  system  complex  frequency  response. 
(AY  =  -  10,  AX  =  -20) 
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CflSErHOMOG   CW  FREQUENCY:    400.0    HZ 

YT:     1000.0    M      XR:-20.0    f1      YR:    990.0    M      ZR:    44.7    h 

DXRi-   1.152    M      DYR:     1.152    M 


Fig.  6.12b  Phase  of  the  overall  system  complex  frequency  response. 
(AY  =  -  10,  AX  =-20) 
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to  the  left  of  the  source.  Remember  also  that  this  geometry  relates  to  the  center  of  the 
receive  array.  Given  this  geometry,  we  expect  that  if  we  move  left  of  center,  which 
corresponds  to  moving  in  a  positive  direction  along  the  X  axis  or  increasing  the  element 
index  m  (see  Fig.  6.4),  the  smaller  the  acoustic  field  would  become.  Conversely,  if  we 
move  to  the  right  of  center,  the  radial  distance  between  the  source  and  the  field  point 
decreases  and  the  acoustic  field  increases.  Also  given  this  geometry,  as  we  go  deeper, 
which  corresponds  to  movement  in  a  positive  direction  along  the  Y  axis  or  increasing  the 
array  element  index  n,  the  radial  distance  between  the  source  and  the  field  point  decreases 
and  the  acoustic  field  increases.  This  behavior  is  clearly  demonstrated  in  Fig.  6.10a. 
Note  that  in  this  explanation,  it  was  assumed  that  the  field  points  did  not  cross  the  dividing 
lines  where  either  AX  or  AY  is  zero. 

Note  the  magnitude  plot  of  Fig.  6.7a,  and  the  way  the  surface  caves  in  by 
approximately  10%  around  the  center  element  along  the  Y  axis,  i.e.,  n  =  0.  This 
corresponds  to  the  center  of  the  receiver  being  at  the  same  depth  as  the  source,  i.e.,  axis- 
axis  propagation.  This  phenomenon  is  caused  by  the  fact  that  our  expression  for  the  ocean 
medium  transfer  function  is  based  on  the  WKB  approximation  and  is  not  exact.  It  was 
observed  in  Chapter  4  that  as  the  propagation  vector  in  the  Y  direction  ky  (y)  approached 
zero,  the  approximate  WKB  solution  is  not  valid  since  the  solution  approaches  infinity. 
This  can  also  be  thought  of  in  terms  of  the  simplifying  assumption  given  by  Eq.  (4.17) 
being  invalid  in  this  region. 

Inspection  of  Eq.  (6.5)  reveals  that  the  integration  with  respect  to  fr  will  cause  the 
ocean  medium  transfer  function  to  be  calculated  for  a  wide  range  of  values  of  the 
propagation  vector  in  the  Y  direction  for  all  spatial  locations.  Why  is  it  only  significant 
when  the  transmitter  and  receiver  are  at  close  to  the  same  depth?  By  thinking  of  the 
problem  from  the  perspective  of  ray  acoustics,  with  the  full-wave  solution  consisting  of  an 
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infinite  summation  of  rays,  it  can  be  seen  that  as  the  transmitter  and  receiver  get  closer 
together  in  depth,  the  eigenray  which  connects  the  two  gets  closer  to  being  horizontal  and 
hence,  closer  to  the  invalid  region  of  the  WKB  approximation.  Since  the  eigenray  and, 
those  close  are  a  major  component  of  the  acoustic  signal  at  the  receiver,  it  is  when  they  get 
close  to  the  invalid  region  that  the  effect  becomes  noticeable. 

This  phenomenon  also  demonstrates  an  important  point  about  ray  acoustics.  In  ray 
acoustics,  rays  other  than  the  eigenrays  or  those  in  very  close  proximity  are  often  ignored 
and  thought  of  as  being  unimportant.  As  a  result,  ray  acoustics  based  on  the  WKB 
approximation  would  be  of  no  use  in  an  along  axis,  isospeed  propagation  problem  since  the 
resulting  acoustic  signal  would  be  infinite.  However,  we  have  seen  that  the  full-wave 
solution  for  an  invalid  region,  by  considering  an  infinite  set  of  rays,  was  still  able  to  give  a 
solution  to  within  10%  of  the  actual  value.  The  full- wave  solution  is  so  called  because  it 
does  not  make  any  assumptions  about  what  is  important,  and  integrates  over  the  entire 
region  of  propagating  and  evanescent  acoustic  wave  contributions.  We  have  shown  that 
this  infinite  integration  allows  the  solution  to  be  equated  to  the  exact  solution  for  the  free- 
space,  isospeed  problem  and  gives  accurate  values  for  the  magnitude  and  phase  of  the 
complex  acoustic  field  at  a  given  field  point. 

G.  VALIDATION  OF  THE  TRANSMIT  ARRAY  FAR-FIELD  DIRECTIVITY 
FUNCTION 

Recall  from  Chapters  2  and  3  the  discussion  regarding  the  far-field  directivity  function, 
or  beam  pattern,  of  a  planar  array  of  point  sources  located  in  the  XY  plane.  If  our 
assumption  about  the  axisymmetric  nature  of  the  far-field  directivity  function  is  correct, 
then  applying  the  analysis  of  the  preceding  sections  in  this  chapter  to  the  overall  system 
complex  frequency  response  given  by  Eq.  (3.14)  yields 
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H(f,r)  =  2jt   I  X  X      cmn(f)HM(f,fr,y0;y)Jo(27tfrrm) 

^0  m  =  -Mf    n  =  -NT 

xe-j2nfY(y0Xy-yn)frdfr;  (6?) 


which  should  be  a  valid  approximate  full-wave  solution  for  our  ocean  medium  propagation 
problem. 

To  validate  Eq.  (6.7),  we  will  compute  the  value  of  the  acoustic  field  at  a  number  of 
field  points  (with  a  RLOS  =  50  m)  due  to  a  3x3  element  transmit  array  and  compare  the 
results  from  exact  calculations  and  the  approximate  full-wave  solution.  The  transmitted 
electrical  signal  will  be  the  Lanczos  smoothed,  rectangular  weighted,  CW  pulse  used 
extensively  throughout  this  chapter.  Recall  that  this  signal  has  a  maximum  frequency 
component  of  640  Hz  and  a  carrier  frequency  of  400  Hz.  For  the  purposes  of  this 
investigation,  we  will  only  consider  the  carrier  frequency  of  the  transmitted  signal. 

First  we  need  to  check  if  we  are  in  the  far  field.  The  intereiement  spacing,  in  both  the 
X  and  Y  directions,  of  the  transmit  array  is  equal  to  one  half  of  the  wavelength 
corresponding  to  the  highest  frequency  component  of  the  transmit  signal.  Ziomek  [Ref. 
3:pp.  1 18-120]  shows  that  this  value  for  the  intereiement  spacing  is  a  necessary  condition 
for  avoiding  grating  lobes  under  all  conditions  of  beam  steering  for  planar  arrays  consisting 
of  MT  x  NT  (odd)  elements.  This  is  the  default  value  used  in  the  computer 
implementation  if  no  user-defined  intereiement  spacing  is  specified,  and  is  the  value  for 
which  we  will  conduct  this  analysis.    Combined  with  Eq.  (2.23),  this  yields  the  condition 


2 

2k  dv-r 

RLOS  >         A1  ,  (6.8) 
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which,  if  satisfied,  places  the  field  point  in  the  far  field  of  the  transmit  array  and  the  far- 
field  directivity  function  is  the  appropriate  description  of  the  beam  pattern  for  this  scenario. 
For  this  case,  the  right-hand  side  of  Eq.  (6.8)  is  equal  to  2.263  m,  which  places  the  RLOS 
of  50  m  well  into  the  far  field 

If  we  further  assume  that  the  beam  pattern  is  not  steered  and  there  is  unit  amplitude 
rectangular  weighting  of  the  array  elements,  it  can  be  shown  that  Eqs.  (2.28)  and  (2.30) 
combine  to  yield 


=  sin[u7t  15/16]   sin[v7t  15/16] 
T^  '  ;      sin[u7c5/16]     sin [vtu  5/16]    '  {     } 


where  the  direction  cosines  u  and  v  are  calculated  from 


and 


AX 
u  =  RLOS  <6'9a> 


AY 
:RLOS-  (610b> 


Values  of  the  magnitude  of  the  overall  system  complex  frequency  response  for  various 
combinations  of  AX  and  AY  are  detailed  in  Table  6.6.  The  exact  value  of  the  acoustic  field 
is  calculated  by  multiplying  the  relevant  value  of  the  far-field  directivity  function,  given  by 
Eq.  (6.8),  with  the  exact  value  of  the  acoustic  field  due  to  a  single  point  source  at  the  same 
RLOS  calculated  from  Eq.  (4.25).  The  simulated  value  of  the  overall  system  complex 
frequency  response,  on  the  other  hand,  is  calculated  using  the  approximate-full  wave 
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solution  given  by  Eq.  (6.7).  As  can  be  easily  observed  from  these  results,  the  value  of  the 
acoustic  field  calculated  using  the  approximate  full-wave  solution  behaves  well  and  is 
within  0.07%  of  the  exact  value  for  a  wide  range  of  direction  cosine  values. 

Recall  from  the  discussion  of  the  results  detailed  in  Table  6.3  that  the  numerical 
integration  of  the  overall  system  complex  frequency  response  was  carried  out  to  sufficient 
terms  to  give  a  maximum  of  a  0.01  %  relative  error.  Note  that  the  errors  detailed  in  Table 
6.6  are  larger  than  this  value.  Does  this  mean  that  there  is  an  error  with  the  assumption 
that  the  beam  pattern  is  axisymmetric  ?  The  answer  is  no.  The  increased  error  is  a  result 
of  the  manner  in  which  Eq.  (6.7)  is  numerically  computed.  To  aid  in  the  convergence  of 
the  numerical  integration  routine,  the  summations  of  Eq.  (6.7)  are  taken  outside  of  the 
integral.  Therefore,  each  integration  within  the  summation,  with  its  associated  error, 
contributes  to  the  total  error  of  the  calculated  overall  system  complex  frequency  response. 

Next,  consider  a  steered  beam  pattern.  The  default  value  for  the  direction  of 
beamsteering  in  the  computer  implementation  is  the  line-of-sight  direction  between  the 
center  of  the  transmit  array  and  the  field  point.  In  other  words,  the  beam  is  steered 
towards  the  field  point.  This  default  value  is  used  if  no  user-defined  direction  is  specified, 
and  this  is  the  case  we  will  investigate  for  this  analysis.  Since  the  beam  is  steered  towards 
the  field  point,  it  can  be  shown  that  the  far-field  directivity  function  given  by  Eq.  (2.28)  in 
the  direction  of  the  field  point  reduces  to  a  constant  equal  to  the  number  of  point  sources  in 
the  array.  In  other  words,  the  array  sources  coherently  add  in  the  direction  of  the  center  of 
the  main  lobe  of  the  beam  pattern.  As  a  result,  the  far-field  directivity  function  of  our 
steered  transmit  array  has  the  value  9  in  the  direction  of  the  field  point.  The  values  of  the 
exact  and  approximate  full-wave  solutions  of  the  acoustic  field  due  to  this  steered  array  are 
detailed  in  Table  6.7.      As  can  be  observed,  like  the  results  contained  in  Table  6.6,  the 
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exact  and  approximate  results  show  close  agreement,  with  errors  being  of  the  same  order  of 
magnitude. 

These  results  confirm  that  our  analysis  relating  to  the  far-field  directivity  function  and 
axisymmetric  behavior  of  a  planar  array  of  point  sources  contained  in  the  previous  chapters 
was  correct.  The  importance  of  this  result  is  that  it  demonstrates  that  the  assumption  made 
in  some  other  ocean  propagation  simulation  programs,  that  axial  symmetry  can  only  be 
used  for  single  point  sources  or  vertical  line  arrays  [Ref.  7]  is  incorrect 

F.   IMPLEMENTATION  OF  NUMERICAL  INTEGRATION 

To  implement  the  integration  of  Eq.  (3.14)  with  respect  to  fr,  the  IMSL10  integration 
routine.  DQDAG,  was  used.  This  is  a  double  precision  numerical  integration  routine 
which  uses  a  globally  adaptive  scheme  based  on  Gauss-Konrad  rules.  Note  that  the 
integrand  of  Eq.(3.14)  contains  both  a  zero-order  Bessel  function  of  the  first  kind, 
implemented  using  the  double  precision  IMSL10  function,  DBSJO,  and  complex 
exponential  terms.  This  means  that  the  integrand  is  oscillatory,  becoming  more  oscillatory 
as  frequency  and  the  range  between  the  source  and  the  field  point  are  increased.  As  a 
consequence,  a  Gauss-Konrad  rule  is  used  with  30-61  points. 

The  oscillatory  nature  of  the  integrand  means  that  it  is  time  consuming  to  compute 
numerically  with  a  accuracy  better  than  0.1%  -  0.5%.  For  our  computer  implementation 
we  desired  a  relative  error  less  than  0.01%,  which  was  not  achievable  by  numerically 
integrating  between  the  appropriate  limits,  regardless  of  how  much  time  was  allowed.  To 
achieve  this  error  criterion,  it  was  necessary  to  subdivide  the  original  regions  of  integration 
into  20  subintervals,  and  to  compute  each  interval  separately  and  sum  the  result. 

It  has  been  mentioned  on  numerous  occasions  that  numerical  integration  is  a  time 
consuming  operation.  To  illustrate  this  fact,  some  typical  run  times  are  detailed  in  Table 
6.8.    Note  that  this  program  has  been  implemented  in  FORTRAN  77  and  was  run  on  the 
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AX 

Dr(f,fx) 

AY 

Dr(f,fY) 

1  H(f,r)  1 
( exact ) 

1  H(f,r)  1 
( simulated ) 

%  error 

0 

3.0000 

25 

2.1111 

1.0079  e-2 

1.0081  e-2 

0.020 

20 

2.4142 

-20 

2.4142 

9.2758  e-3 

9.2795  e-3 

0.040 

10 

2.8478 

10 

2.8478 

1.2907  e-2 

1.2903  e-2 

-  0.031 

-10 

2.8478 

-30 

1.7654 

8.0015  e-3 

8.0038  e-3 

0.029 

35 

1.3902 

10 

2.8478 

6.3008  e-3 

6.3050  e-3 

0.067 

Table  6.6   Acoustic  field  resulting  from  an 
unsteered,  3x3  transmit  array. 


AX 

AY 

1  H(f,r)  1 
( exact ) 

1  H(f,r)  1 
( simulated ) 

%  error 

10 
20 

10 
-20 

1.4324  e-2 
1.4324  e-2 

1.4318  e-2 
1.4320  e-2 

-  0.042 

-  0.028 

Table  6.7  Acoustic  field  resulting  from  a  3  x  3  transmit  array 
steered  in  the  direction  of  the  field  point. 
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Naval  Postgraduate  School  IBM  3033/4381  Network.  It  should  be  noted  that  the  run 
times  shown  in  Table  6.8  are  for  a  frequency  of  400  Hz.  Note  that  the  times  detailed  in 
Table  6.8  are  for  the  simplest  case  possible.  The  phase  integral  term  of  the  ocean  medium 
transfer  function  has  a  closed-form  expression  for  the  isospeed  medium  case,  which  means 
that  no  numerical  integration  or  approximation  of  the  phase  integral  has  to  be  made. 
Hence,  computation  time  for  this  component  of  the  integrand  is  at  a  minimum.  The 
computation  time  becomes  very  significant  when  dealing  with  receive  and  /  or  transmit 
arrays  rather  than  single  elements.  This  is  the  main  reason  that  most  of  our  analysis  has 
been  carried  out  for  small  values  of  RLOS. 


RLOS 

Computation 

(m) 

Time  ( min:sec ) 

10 

2:01 

100 

2:20 

1000 

3:23 

10000 

12:18 

Table  6.8  Simulation  run  times  for  a  single  point  source 
to  a  single  point  receiver  problem  at  400  Hz. 
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VII.   COMPUTER  SIMULATION  RESULTS  FOR  THE 
SINGLE  GRADIENT  MEDIUM  CASE 

A.    EVALUATION  OF  THE  WKB  PHASE  INTEGRAL 

The  results  of  the  previous  chapter  validated  the  development  carried  out  in  Chapters  2 
through  4  for  the  case  of  a  free-space,  isospeed  medium.  The  next  step  in  the  validation 
process  is  to  apply  the  analysis  to  a  more  complex  scenario,  i.e.,  the  free-space,  single 
gradient  case.  For  the  purposes  of  this  thesis  we  will  characterize  the  sound-speed  profile 
as  being  linear,  with  the  speed  of  sound  being  given  by 

c(y)  =  cYref+  g(  y  -  YrefX  (7-1) 

where  g  is  the  sound-speed  gradient,  and  cYref  *s  tne  speed  of  sound  at  the  depth  yRgp 
The  linear  approximation  of  the  sound-speed  profile  was  selected  as  it  allows  for  the  easy 
description  of  a  more  complex  medium  by  means  of  linear  segments.  The  aim  of  this 
chapter  is  to  validate  the  analysis  for  a  single  linear  segment  before  the  computer 
implementation  is  applied  to  a  multiple  linear  segment  sound-speed  profile. 

Since  the  preceding  development  of  the  overall  system  complex  frequency  response 
given  by  Eq.  (3.14)  assumed  a  free-space  problem,  all  we  need  to  do  is  apply  the  condition 
of  a  single,  linear,  sound-speed  profile  to  the  ocean  medium  transfer  function,  given  by 


HM(f,fr,yo:y)=  -Ja-ei'^^WoXy-yJ        (4.34) 

2kY(yo)V|kY(y)| 
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where 


and 


2     2   1/2        2  2 

kY(y)  =  ±27c[(f/c(y))   -f{]      ,    i7<(f/c(y)),  (7.2a) 


2  l}11        2  2 

kY(y)  =  +  j27u[fr  -(f/c(y))1      ,    fr>(f/c(y))z  (7.2b) 


For  the  form  of  the  sound-speed  profile  given  by  Eq.  (7.1),  we  were  unable  to  find  a 
closed  form  solution  to  the  phase  integral  in  Eq.  (4.34).  This  is  in  marked  contrast  to  the 
results  of  the  previous  chapter,  where  the  phase  integral  had  a  closed  form  expression 
which  resulted  in  the  ocean  medium  transfer  function  for  the  isospeed  medium  requiring  no 
numerical  integration  in  its  evaluation. 

The  method  required  to  evaluate  the  phase  integral  of  Eq.  (4.34)  is  of  critical 
importance.  Recall  from  Chapter  3  that  the  expression  for  the  overall  system  complex 
frequency  response,  given  by  Eq.  (3.14),  involved  the  evaluation  of  an  integral  whose 
integrand  contained  the  ocean  medium  transfer  function.  In  other  words,  there  exists  an 
integration  within  an  integration.  If  a  numerical  integration  routine,  such  as  the  IMSL10 
routine,  DQDAG,  discussed  in  the  previous  chapter,  is  used  to  evaluate  this  phase  integral, 
then  the  resulting  computational  time  involved  is  unmanageable.  As  an  example,  the 
computer  simulation  was  run  for  a  single  frequency  component  at  a  RLOS  of  100  m,  and  a 
single  point  source  and  a  single  point  receiver.  It  was  found  for  this  case  that  the 
maximum  allowable  processing  time  of  one  hour  on  the  Naval  Postgraduate  School  IBM 
3033/4381  Network  was  attained  before  a  solution  could  be  computed.  Consequently,  we 
required  an  alternative  method  to  numerical  integration  to  evaluate  the  integral. 
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As  an  approximate  solution  to  the  phase  integral,  we  expressed  the  propagation  vector 
components  in  the  Y  (depth)  direction,  given  by  Eqs.  (7.2a)  and  (7.2b),  in  terms  of  their 
binomial  series.  The  reason  for  doing  this  was  that  the  integral  of  each  of  the  terms  of  the 
binomial  series  expansion,  for  both  forms  of  the  propagation  vector  component  ky  (y),  has 
a  closed  form  expression.  As  a  result,  the  phase  integral  of  the  ocean  medium  transfer 
function  can  be  replaced  by  a  converging  series.  Taking  the  integral  of  the  binomial  series 
expansion  of  Eqs.  (7.2a)  and  (7.2b),  it  can  be  shown  that 

(  .     ,  .  .       ,  2rcf  r  .  ,  ,  N ,      1  /  frc(y)  \   l      l  /  frc(y)  \    i 
J    kY(y)dy  =  ±—  |_ln[c(y)]-2(-T-)2-8(-f-)    4 

/frc(y)\6i       5    /frc(y)\8i  1 

(-J")  6  "128  ("J-)    8  J  (7-3> 


16  ' 

for  the  propagating  waves  described  by  Eq.  (7.2a),  and 

3 


/kY(y)dy  =  ,j2lt{-fry-l['(^)   +Kriy)) 

+  4(f^y))   J  +  T§8(f^yy)    T+        ]    }     <74) 

for  the  evanescent  waves  described  by  Eq.  (7.2b).     The  number  of  terms  used  in  Eqs. 
(7.3)  and  (7.4)  depends  on  the  desired  accuracy  for  the  phase  integral.     Note  that  this 

method  of  determining  the  phase  integral  is  inefficient  for  small  values  of  ky  (y).    This  is 

2  2 

due  to  the  fact  that  for  small  values  of  ky  (y),  fr  =  (f  /c(y))  .    As  a  result,  the  binomial 

expansion,  and  hence,  the  phase  integral  converge  very  slowly. 
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Recall  from  Chapter  3  that  the  plus  (minus)  sign  in  Eq.  (7.3)  is  chosen  whenever 
y  >  yG  (y  <  yo)  and  corresponds  to  the  field  point  being  deeper  (shallower)  than  the  source 
point,  i.e.,  downward  (upward)  traveling  waves.  The  minus  (plus)  sign  in  Eq.  (7.4) 
corresponds  to  the  plus  (minus)  sign  of  Eq.  (7.3)  in  order  to  generate  evanescent  waves. 
Recall  that  the  full-wave  solution  can  be  thought  of  as  an  infinite  summation  of  rays.  If  we 
think  of  rays  from  a  single  omnidirectional  point  source  propagating  in  a  medium 
characterized  by  a  single,  linear,  sound-speed  gradient,  it  can  easily  be  visualized  that  some 
rays  will  never  reach  the  receiver  depth  while  others  will  reach  the  receiver  depth  twice. 
This  situation  is  shown  in  Fig.  7.1  where  it  has  been  assumed  that  the  receiver  is  below  the 
source  and  the  sound-speed  gradient  is  positive. 

This  is  an  important  concept  when  we  consider  the  way  in  which  the  overall  system 
complex  frequency  response  is  calculated.  Since  the  overall  system  complex  frequency 
response  is  an  integral  with  respect  to  fr  and  is  a  function  of  frequency,  the  WKB  phase 
integral  contained  in  the  integrand  is  essentially  the  phase  change  of  a  particular  ray, 
defined  by  the  values  of  f  and  fr,  traveling  between  the  source  and  the  receiver  depth.  For 
each  of  these  rays,  we  need  to  decide  how  it  behaves  and  determine  a  strategy  for 
determining  the  value  of  the  phase  integral. 

This  ray  path  problem  was  not  encountered  when  dealing  with  the  isospeed  medium, 
since  for  this  case  the  rays  travel  in  straight  lines  as  shown  in  Fig.  7.2.  The  situation 
shown  in  Fig.  7.2  assumes  that  the  receiver  is  below  the  source  and,  as  can  be  seen,  all  the 
downward  traveling  waves  reach  the  receiver  depth  once  in  an  isospeed  medium.  As  a 
consequence,  only  the  closed  form  expression  of  the  phase  integral  for  downward  traveling 
waves  needs  to  be  used  in  the  evaluation  of  the  ocean  medium  transfer  function  for 
propagating  waves.  Similarly,  only  one  form  of  the  phase  integral  expression  needs  to  be 
used  for  the  evanescent  waves. 
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( cross  range ) 


Fig.  7.1  Ray  paths  in  a  medium  characterized  by 
a  single,  positive,  sound-speed  gradient. 


( cross  range ) 


Fig  7.2  Ray  paths  in  an  isospeed  medium. 
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Inspection  of  Fig.  7.1  shows  the  various  groups  of  rays  that  will  be  encountered. 
Rays  A  and  B  are  handled  in  the  same  manner.  Since  neither  of  these  rays  is  refracted 
back  to  the  receiver  depth  before  reaching  the  horizontal  range  of  the  receiver,  the  phase 
integral  is  calculated  for  the  downward  traveling  portion  of  the  rays.  For  rays  A  and  B, 
the  integral  is  calculated  from  the  source  depth  y0  meters  to  the  first  attainment  of  the 
receiver  depth  y  meters,  i.e.,  from  y0  to  points  1  and  2,  respectively. 

Ray  C,  unlike  rays  A  and  B,  is  refracted  back  to  the  receiver  depth  before  the  horizontal 
range  of  the  receiver.  As  is  illustrated  in  Fig.  7. 1 ,  the  ray  has  both  downward  and  upward 
portions  that  need  to  be  taken  into  consideration.  Consequently,  the  phase  integral  is 
calculated  for  a  downward  traveling  ray  from  the  source  to  the  turning  point  at  point  3,  and 
for  an  upward  traveling  ray  from  the  turning  point  back  to  the  receiver  depth  at  point  4. 
This  essentially  means  that  the  phase  integral  has  been  subdivided  into  two  separate 
segments,  dependent  on  which  form  of  the  propagation  vector  given  by  Eq.  (7.2a)  is  valid. 
Brekhovkikh  [Ref.  6:pp.  66-67]  also  shows  that  such  a  ray  encounters  a  the  phase  shift  of 
-7t/2  radians  at  the  turning  depth,  i.e.,  point  3.  After  passing  through  the  turning  point,  the 
upward  traveling  ray  keeps  the  -n/2  phase  shift. 

The  final  ray  type  to  be  considered  is  that  characterized  by  ray  D.  As  can  be  seen  from 
Fig.  7.1,  this  propagating  ray  never  reaches  the  receiver  depth  since  it  reaches  a  turning 
point  above  the  receiver  depth  at  point  5,  and  becomes  an  upward  traveling  ray.  To 
evaluate  the  phase  integral,  we  first  consider  the  downward  traveling  portion  of  the  ray 
down  to  the  depth  of  the  turning  point.  We  then  need  to  evaluate  the  integral  from  this 
point  down  to  the  receiver  depth.  There  is  clearly  no  propagating  wave  component  going 
down  to  the  depth  of  the  receiver  from  the  turning  point,  so  any  acoustic  field  reaching  the 
depth  of  the  receiver  from  this  point  must  be  evanescent.  Accordingly,  to  evaluate  the 
phase  integral  down  to  the  receiver  depth  we  consider  the  integrand  to  be  the  complex 
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propagation  vector  component  given  by  Eq.  (7.4).  Recalling  the  discussion  regarding  the 
behavior  of  ray  C  at  a  turning  point,  this  evanescent  component  can  be  thought  of  as  what 
gives  rise  to  the  -rc/2  phase  shift  in  the  propagating  portion  of  the  ray  reflected  at  the  turning 
point.  In  other  words,  when  the  ray  goes  through  a  turning  point,  it  gives  rise  to  a  phase 
shifted,  reflected  propagating  component  and  a  evanescent  component. 

Recall  from  the  previous  chapter  the  importance  of  evanescent  waves.  As  was  done 
for  the  isospeed  medium  case,  we  will  use 

-1.2(f/c0) 
H(f,r)  =  2k   I  HM  (f,fr,y0;y)  JoG^m )  e"j27tfY(yoXy  "  yn)fr  dfr.    (6.5) 

Jo 

as  our  approximate  full-wave  solution  to  the  ocean  propagation  problem  which  means  that 
we  will  be  integrating  with  respect  to  fr  up  to  a  limit  of  1.2  (f  /  cG).  Since  this  integration 
limit  gave  very  good  results  for  the  isospeed  medium  case,  it  is  assumed  that  it  will  also  be 
applicable  to  the  new  single,  linear  gradient  problem  since  the  sound-speedprofile  is  slowly 
changing  and  does  not  exhibit  any  erratic  behavior.    Note  that  the  region  of  integration 

(f/c0)<fr<1.2(f/c0)  (7.4) 

describes  evanescent  waves  resulting  from  the  source  condition  and  should  not  be  confused 
with  an  evanescent  wave  that  results  from  a  ray  being  refracted  through  a  turning  point. 

B.  RESULTS 

As  discussed  in  the  preceding  section,  the  accuracy  of  the  binomial  series 
approximation  depends  on  the  number  of  terms  used  and  for  small  values  of  kY(y)  the 
number  of  terms  required  becomes  very  large.    Just  as  for  the  situation  in  which  we  tried 
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to  evaluate  the  phase  integral  using  a  numerical  integration  routine,  if  we  specify  a 
requirement  that  the  phase  integral  be  accurate  to  within  0.01%  and  evaluate  the  binomial 
expansion  approximation  for  all  the  values  of  fr  used  to  numerically  integrate  Eq.  (6.5),  the 
computation  time  becomes  unmanageable.  In  this  case  it  was  again  found  that  the 
maximum  allowable  processing  time  of  one  hour  on  the  Naval  Postgraduate  School  IBM 
3033/4381  Network  was  attained  before  a  solution  could  be  computed. 

As  a  consequence,  it  was  decided  to  modify  the  region  of  integration  in  the  evaluation 
of  Eq.  (6.5).  Instead  of  integrating  over  the  two  regions  0  <  fr  <  (f /c0)  and 
(f  /c0)  <  fr  <  1.2  (f /c0),  we  calculated  the  approximate  full-wave  solution  of  Eq.  (6.5) 
using  the  following  modified  regions:  0  <  f  r  <  LIM  (f  /  c0 )  for  initially  propagating  waves 
and  (1/LIM)  (f /c0)  <  fr  <  1.2  (f/c0)  for  initially  evanescent  waves.  In  these 
expressions,  LIM  is  a  constant  less  than  unity  that  basically  allows  us  to  trade  computation 
time  against  accuracy.  After  some  investigation  of  the  rate  of  convergence  of  the  binomial 
series  approximation  of  the  phase  integral,  it  was  found  that  solutions  could  be  obtained  for 
a  scenario  in  which  the  RLOS  was  100  m,  the  transmitter  was  a  single  point  source,  and 
the  receiver  was  a  3  x  3  array  of  point  receivers,  if  the  value  of  LIM  did  not  exceed  0.995. 
It  can  easily  be  seen  that  the  closer  the  value  of  LIM  is  to  unity,  the  better  our  approximate 
full-wave  solution  should  be.  Accordingly,  we  conducted  simulations  for  LIM  equal  to 
0.99  and  0.995. 

From  our  observations  in  the  previous  chapter,  there  are  two  forms  of  output  that  are 
useful  in  determining  the  performance  of  our  full-wave  approximation.  First,  there  is  the 
tabulation  of  the  magnitude  and  phase  of  the  complex  acoustic  field  at  the  locations  of  the 
receive  array  elements.  Recall  that  the  usefulness  of  this  form  of  output  was  that  it  allows 
for  easy  verification  that  the  magnitude  of  the  complex  acoustic  field  is  behaving  correctly. 
For  a  given  source  /  receiver  geometry,  it  is  an  easy  matter  to  determine  if  the  the  magnitude 
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of  the  complex  acoustic  field  is  increasing  or  decreasing  according  to  theory  as  we  go  from 
one  element  to  the  next.  Second,  there  is  the  output  from  the  adaptive  beamforming 
algorithm,  used  in  the  previous  chapter,  processing  the  time-domain  electrical  signal  at  each 
element  of  the  receive  array.  This  algorithm  primarily  uses  the  phase  of  the  received  signal 
in  its  processing.  Therefore,  it  validates  how  accurate  the  phase  of  the  predicted  complex 
acoustic  field  is.  Consequently,  analysis  of  both  these  outputs  provides  a  check  of  the 
accuracy  of  both  the  amplitude  and  phase  of  the  complex  acoustic  field  at  a  specified  field 
point. 

Computer  simulation  results  were  generated  for  two  values  of  the  constant  LIM,  as 
previously  discussed,  in  conjunction  with  a  source  /  receiver  geometry  of  the  form  shown 
in  Fig.  7.1.  The  source  was  placed  at  a  depth  of  1,000  m,  the  receiver  was  below  at  a 
depth  of  1,050  m,  and  the  RLOS  was  100  m.  The  sound-speed  profile  of  the  ocean 
medium  was  characterized  by  the  linear  sound-speed  profile  given  by  Eq.  (7.1),  where 
VREF=  1>000  m,  Cyref=  1»475  m/sec  and  the  gradient  had  the  value  0.017  1/sec.  The 
aim  of  investigating  these  particular  values  of  LIM  was  to  determine  whether  the  small 
region  of  integration  that  was  discarded  was  significant. 

The  values  of  the  magnitude  and  phase  (in  degrees)  of  the  acoustic  field  for  four  test 
cases  are  shown  in  Tables  7.1  through  7.4.  Note  that  the  first  value  is  the  magnitude  and 
the  second  value  is  the  phase.  Given  the  placement  of  the  receive  array  relative  to  the 
source,  i.e.,  XR  =  0  and  YR  =  1050  m,  which  is  below  the  source,  we  would  expect  the 
magnitude  of  the  acoustic  field  to  decrease  with  increasing  n.  This  corresponds  to  reading 
down  the  table  or  equivalently  increasing  the  depth  of  the  receive  element  or  field  point. 
We  would  also  expect  the  value  of  the  acoustic  field  to  decrease  with  increasing  absolute 
value  of  m.  This  corresponds  to  reading  across  the  table  out  from  the  center  or 
equivalently  increasing  the  displacement  in  the  X  direction  of  the  element  or  field  point. 
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We  expect  the  magnitude  of  the  acoustic  field  to  decrease  in  both  of  these  situations  since 
their  effect  is  to  increase  the  radial  range  between  the  source  and  the  receive  element ,  or 
field  point  Note  that  these  are  the  same  results  as  expected  for  the  isospeed  medium  case. 
This  is  because  the  RLOS  is  a  small  value,  and  hence,  the  propagating  rays  have  not  had 
much  of  an  opportunity  to  refract  Therefore,  a  crude  approximation  is  that  the  medium  is 
homogeneous  over  this  short  range. 

Case  A,  shown  in  Table  7.1,  evaluates  the  overall  system  complex  frequency  response 
over  the  modified  region  of  the  initially  propagating  waves.  In  other  words,  we  only 
integrate  over  the  region  0  <  f  r  <  LQM  (f  /  c0 )  and  ignore  the  region  of  initially  evanescent 
waves.  Case  B.  on  the  other  hand,  shown  in  Table  7.2,  is  the  same  as  case  A  except  that 
it  also  includes  the  modified  region  of  intitially  evanescent  waves.  Both  cases  are 
computed  for  LIM  =  0.99.  It  can  be  easily  observed  from  these  two  tables  that  the 
magnitude  of  the  acoustic  field  does  not  behave  exactly  as  expected.  In  both  cases,  the 
magnitude  of  the  acoustic  field  decrease  as  we  read  down  the  table  as  expected.  However, 
in  both  cases,  the  second  and  third  rows  have  the  magnitude  increasing  as  we  move  away 
from  m  =  0. 

Note  that  the  magnitude  and  phase  of  the  complex  acoustic  field  is  very  similar  for  both 
these  cases  despite  the  inclusion  of  evanescent  waves  in  Case  B.  At  first,  this  may  appear 
to  contradict  the  results  of  the  previous  chapter.  For  the  isospeed  medium  the  inclusion  of 
evanescent  waves  had  a  dramatic  effect  on  the  value  of  the  complex  acoustic  field.  So 
what  causes  this  apparently  contradictory  behavior?  The  difference  between  the  observed 
results  for  the  isospeed  and  single  sound-speed  gradient  cases  is  due  to  the  region  of 
evanescent  waves  considered.  The  evanescent  waves  in  the  region  very  close  to 
fr  =  (f /c0)  are  by  far  the  most  dominant  ones.  Accordingly,  the  modified  region  of 
evanescent  waves  ignores  the  most  dominant  ones.    So  why  the  upper  limit  of 
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n/m 

1 

0 

-1 

-1 

8.17543  e-4 
-  165.452 

8.18198  e-4 
-  164.725 

8.17543  e-4 
-  165.452 

0 

8.12795  e-4 
139.886 

8.12545  e-4 
140.616 

8.12795  e-4 
139.886 

1 

7.98127  e-4 
83.940 

7.97160  e-4 
84.627 

7.98127  e-4 
83.940 

Table  7.1   Case  A  :   LIM  =  0.99,  initially  evanescent  waves  not  included. 


n/m 

1 

0 

-1 

-1 

8.17524  e-4 
-  165.451 

8.18179  e-4 
-  164.725 

8.17524  e-4 
-  165.451 

0 

8.12783  e-4 
139.885 

8.12534  e-4 
140.615 

8.12783  e-4 
139.885 

1 

7.98128  e-4 
83.939 

7.97161  e-4 
84.626 

7.98128  e-4 
83.939 

Table  7.2   Case  B  :   LIM  =  0.99,  initially  evanescent  waves  included. 
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fr  =  1.2  (f/c0)?  This  upper  limit  becomes  significant  when  considering  very  short 
ranges,  i.e..  ranges  less  than  100  m.  At  these  short  ranges  the  region  of  significant 
evanescent  waves  increases  rapidly  with  decreasing  range.    The  specific  upper  limit  of 

fr  =  1.2  (f /c0 )  was  chosen  to  give  accurate  values  for  the  acoustic  field  calculated  from 
the  approximate  full-wave  solution  down  to  a  range  of  10  m. 

If  our  reasoning  has  been  correct,  then  we  would  expect  our  values  for  the  computed 
complex  acoustic  field  to  improve  as  we  increased  the  value  of  LIM  closer  to  unity.  Case 
C,  which  ignores  the  modified  region  of  evanescent  waves,  and  Case  D.  which  includes 
the  modified  region  of  evanescent  waves,  were  both  computed  for  LIM  =  0.995  and  are 
shown  in  Tables  7.3  and  7.4.  respectively  As  can  be  seen  from  the  values  in  these  tables. 
the  same  phenomena  is  demonstrated  as  was  evident  when  comparing  Tables  7.1  and  7.2. 
In  fact,  it  can  be  easily  seen  that  in  a  §  had  minimal 

effect  on  the  value  of  the  complex  acoustic  field.  Does  this  mean  that  our  explanation  of 
the  importance  of  the  evanescent  waves  in  the  region  very  close  to  fr  =  (f/c0)  is 
incorre.  :er  the  isosrx:  _  .::.  .    ipterandscc 

happens  when  we  apply  these  modified  regions  of  integration. 

Consider  the  same  source  /  receiver  geometry  as  we  have  been  using  for  this 
investigation,  but  apply  the  approximate  full-wave  solution  to  an  isospeed  problem,  where 
the  constant  speed  of  sound  c0  is  1.475  vafst  :he  case  in  the  previous  chapter. 

Applying  the  modified  regions  of  integration,  with  LIM  =  0.999,  to  the  ocean  medium 
transfer  function  given  by  Eq.  (6.5  j  results  in  the  values  of  the  complex  acoustic  field 
shown  in  Table  7.5.  Note  the  way  that  the  magnitudes  of  the  acoustic  field  do  not  behave 
as  predicted.  In  the  third  row  the  magnitude  of  the  acoustic  field  increases  as  we  move 
away  from  m  =  0.  Also,  note  that  the  magnitudes  of  the  acoustic  field  increase  with 
increasing  depth.    Both  of  these  results  are  the  opposite  from  what  we  expect.    Comparing 
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n/m 

1 

0 

-1 

-1 

8.29299  e-4 
-  162.942 

8.28427  e-4 
-  162.348 

8.29299  e-4 
-  162.942 

0 

7.94735  e-4 
142.244 

7.94784  e-4 
142.804 

7.94735  e-4 
142.244 

1 

7.61585  e-4 
84.539 

7.62560  e-4 
85.126 

7.61585  e-4 
84.539 

Table  7.3   Case  C  :   LIM  =  0.995,  initially  evanescent  waves  not  included. 


n/m 

1 

0 

-1 

-1 

8.291 10  e-4 
-  162.938 

8.28241  e-4 
-  162.344 

8.291 10  e-4 
-  162.938 

0 

7.94604  e-4 
142.237 

7.94653  e-4 
142.797 

7.94604  e-4 

142.237 

1 

7.61599  e-4 
84.529 

7.62572  e-4 
85.116 

7.61599  e-4 
84.529 

Table  7.4  Case  D  :   LIM  =  0.995,  initially  evanescent  waves  included. 
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Tables  7.5  and  6.5,  we  observe  that  there  is  a  significant  difference  between  the  two  sets  of 
results,  yet  the  only  difference  between  the  two  was  that  a  modified  region  of  integration 
was  used  to  calculate  the  values  in  Table  7.5.  Recall  that  the  value  of  LIM  for  the  values 
calculated  in  Table  7.5  was  0.999,  that  is,  very  close  to  unity.  These  results  reaffirm  our 
earlier  conclusions  that  the  evanescent  waves  in  the  region  very  close  to  fr  =  (f /c0) 
contribute  very  significantly  to  the  overall  value  of  the  integral.  This  implies  that  around 
this  value  of  fr,  there  may  be  a  region  of  relatively  constant  phase  in  the  integrand  of  the 
overall  system  complex  frequency  response. 

The  outputs  from  the  LMS  adaptive  beamforming  algorithm  for  Cases  A  through  D  are 
shown  in  Table  7.6.  Recall  the  geometry  illustrated  in  Fig.  6.3.  The  angle  definitions  are 
still  valid  with  the  exception  that  they  relate  to  the  local  angles  of  arrival  of  the  wave  at  the 
receive  transducer.  In  other  words,  for  the  isospeed  medium  case  the  local  angle  of  arrival 
was  also  the  direction  of  the  source  since  all  rays  traveled  in  straight  lines.  From  the 
values  for  all  four  cases  we  see  that  the  range  estimates  RLOS  (X)  are  in  error  by  an  order 
of  10  -  15  %  while  the  estimates  RLOS  (Y)  are  in  error  by  up  to  65  %.  This  phenomena 
was  also  observed  for  the  isospeed  case  when  the  effects  of  evanescent  were  ignored.  For 
that  case,  it  was  argued  that  since  the  overall  system  complex  frequency  response  is 
axisymmetric,  then  the  relative  behavior  of  the  acoustic  field  in  the  X  direction  should  be 
minimally  affected  by  the  behavior  of  the  acoustic  field  in  the  Y  direction.  Consequently, 
if  there  is  a  problem  with  the  description  of  the  full-wave  solution  in  the  Y  direction,  as  in 
this  case,  the  value  of  RLOS  (X)  should  still  be  fairly  accurate,  which  is  the  result  that  we 
observed. 

The  angular  outputs  shown  in  Table  7.6  agree  very  closely  with  theory.  Even  though 
the  angles  are  the  predicted  local  angle  of  arrival,  for  our  case  where  the  RLOS  is  small 
(i.e.,  100  m),  the  effects  of  refraction  will  have  been  small  and  the  local  angle  of 
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n/m 

1 

0 

-1 

-1 

7.33866  e-4 
-  169.589 

7.34792  e-4 
-  168.589 

7.33866  e-4 
-  169.321 

0 

7.75856  e-4 
131.862 

7.77764  e-4 
132.604 

7.75856  e-4 
131.862 

1 

8.30442  e-4 
76.103 

8.29309  e-4 
76.781 

8.30442  e-4 
76.103 

Table  7.5  Case  E  :   LIM  =  0.999,  initially  evanescent  included. 


CASE 

LEVI 

RLOS  (X) 
(m) 

est 

RLOS  (Y) 
(m) 

est 

0 

(deg) 
est 

(deg) 
est 

A 

0.99 

90.690 

75.577 

29.471 

270.00 

B 

0.99 

90.688 

75.628 

29.471 

270.00 

C 

0.995 

115.748 

34.343 

30.015 

270.00 

D 

0.995 

115.720 

34.438 

30.019 

270.00 

Table  7.6  Output  from  adaptive  LMS  beamforming  algorithm. 
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arrival  should  be  very  close  to  the  angles  describing  the  position  of  the  source  relative  to  the 
receiver.  Any  effect  due  to  refraction  for  the  given  geometry  and  sound-speed  profile  will 
be  to  decrease  the  angle  0  slightly  and  leave  *F  unaffected.  However,  our  results  are  too 
inconsistent  to  tell  if  there  is  a  noticeable  decrease  in  the  value  of  0,  but  they  do 
demonstrate  that  *F  remains  unchanged. 

The  results  obtained  from  this  investigation  of  a  medium  characterized  by  a  linear 
sound- speed  profile  with  a  single  gradient  have  essentially  been  inconclusive.  It  was 
found  that  the  evaluation  of  the  phase  integral  in  the  ocean  medium  transfer  function  based 
on  the  WKB  approximation  is  of  crucial  importance  to  the  performance  of  the  approximate 
full-wave  solution.  If  the  linear  systems  approach,  and  the  coupling  equations  developed 
in  this  thesis  are  to  be  successfully  applied  to  the  WKB  approximation  of  the  ocean 
medium,  we  require  a  fast,  accurate  method  of  determining  the  value  of  the  phase  integral 
to  be  able  to  accurately  compute  the  value  of  the  complex  acoustic  field  at  some  given 
spatial  location. 

The  method  of  representing  the  phase  integral  in  terms  of  a  power  series  expansion, 
while  better  that  using  numerical  integration  techniques,  was  inappropriate  since  there  is  a 
major  contribution  to  the  value  of  the  acoustic  field  in  the  region  where  this  expansion 
converges  very  slowly.  The  use  of  the  series  expansion  allowed  approximate  values  of  the 
acoustic  field  to  be  calculated,  whereas  the  use  of  the  IMSL10  numerical  integration  routine 
DQDAG  took  too  much  computing  time  and  no  useful  results  could  be  obtained. 
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Vm.  CONCLUSIONS  AND  RECOMMENDATIONS 

It  appears  from  our  preliminary  results  that  our  approach  to  modelling  the  ocean 
medium  pulse  propagation  problem,  which  incorporates  both  linear  systems  theory  and  the 
physics  of  acoustic  propagation,  is  a  valid  one. 

A  signal  generation  scheme  based  on  a  truncated  Fourier  series  and  complex  envelope 
representation  of  narrowband  signals  was  implemented  and  shown  to  work  correctly.  The 
scheme  is  capable  of  generating  CW  and  LFM  pulses  with  rectangular,  Hamming  or 
Hanning  amplitude  weighting  suitable  for  conversion  into  a  transmitted  acoustic  signal. 

We  have  shown  that  the  beam  pattern  of  the  transmit  array  does  not  violate  axial 
symmetry.  This  is  a  result  of  the  theory  of  linear  superposition,  which  allows  us  in  the 
case  of  our  linear  system  model,  to  decompose  the  beam  pattern  into  a  sum  of 
omnidirectional  point  sources. 

The  model  has  been  extensively  tested  and  validated  for  the  case  of  a  free-space 
isospeed  medium.  It  was  found  that  the  formal  solution  to  the  propagation  problem  could 
be  modified  to  include  only  a  small  portion  of  the  region  of  evanescent  waves.  This 
approximate  full-wave  solution  has  demonstrated  that  it  gives  extremely  accurate  results  for 
radial  ranges  between  source  and  receiver  down  to  10  m.  It  was  also  shown  that  the 
formal  solution  was  inaccurate  in  the  situation  where  the  source  and  field  point  were  at  or 
very  near  the  same  depth.  It  was  concluded  that  this  was  due  to  the  nature  of  the  WKB 
approximation  rather  than  a  problem  with  the  system  model.  It  is  recommended  that  a 
modified  WKB  approximation  be  used  in  the  vicinity  of  a  turning  point,  which  involves 
Airy  functions.  This  modified  solution  should  be  investigated  in  an  attempt  to  determine  if 
it  improves  the  model  performance  for  this  source  /  receiver  geometry. 
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Inconclusive  results  were  obtained  about  the  model's  behavior  for  an  ocean  medium 
characterized  by  a  linear  sound- speed  profile  with  a  single  gradient.  This  was  due  to  the 
difficulties  encountered  in  evaluating  the  phase  integral  term  of  the  ocean  medium  transfer 
function  based  on  the  WKB  approximation.  For  this  case  we  were  unable  to  determine  a 
closed  form  expression  for  the  phase  integral.  Numerical  integration  techniques  were  too 
slow  and  our  derived  approximation  to  the  integral  was  unsuited  to  the  numerical  behavior 
of  the  integrand  and  its  effect  on  the  behavior  of  the  approximate  full- wave  solution  to  the 
propagation  problem.  Further  investigation  needs  to  be  carried  out  on  how  to  quickly  and 
accurately  evaluate  this  integral.  Some  recommended  avenues  for  further  investigation 
are : 


to  determine  whether  a  different  description  of  the  sound-speed  profile  will  yield  a 
closed  form  expression  for  the  phase  integral, 

to  find  a  better  approximation  of  the  phase  integral  than  the  binomial  series 
expansion. 


There  is  also  a  question  regarding  the  strategy  required  to  implement  the  phase  integral. 
From  a  ray  acoustics  viewpoint,  this  integral  represents  the  phase  change  experienced  by  a 
ray  as  it  travels  from  the  source  to  the  receiver  depth.  For  the  inhomogeneous  medium, 
these  rays  will  be  refracted  and  there  exists  the  situation  were  the  ray  doesn't  reach  the 
receiver  depth  at  all,  or  reaches  it  more  than  once.  How  should  the  integral  be  evaluated  in 
these  circumstances?  Because  of  the  inconclusive  nature  of  our  computer  simulation 
results,  we  were  unable  to  validate  or  disprove  our  scheme  to  compute  the  phase  integral. 

It  was  demonstrated  that  the  outputs  from  the  computer  simulation,  i.e.,  the  magnitude 
and  phase  of  the  complex  acoustic  field  as  a  function  of  frequency  and  spatial  coordinates, 
and  the  time-domain  output  electrical  signal  from  each  element  in  a  receive  array  were 
implemented  correctly.     While  no  matched-field  processing  was  done  using  this  first 
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output,  it  was  demonstrated  that  the  plotting  of  the  magnitude  of  the  complex  acoustic  field 
was  a  powerful  aid  in  helping  to  visualize  what  the  acoustic  field  looked  like.  The  time- 
domain  signals  were  processed  by  an  external  frequency  domain  adaptive  beamforming 
program  and  the  data  obtained  was  very  useful  in  validating  the  performance  of  the  model. 
This  also  demonstrated  that  generated  data  could  be  successfully  passed  to  an  external  file, 
to  be  processed  by  an  ancillary  program  of  our  choosing.  Thus,  the  generated  data  is 
independent  of  whatever  receive  array  signal  processing  we  wish  to  carry  out 

It  has  been  shown  that  an  approximate  solution  to  the  linear  wave  equation  can  be 
implemented  as  an  ocean  medium  transfer  function.  However,  except  when  closed  form 
solutions  or  short  recursive  computation  schemes  exist  for  all  terms  in  the  ocean  medium 
transfer  function,  we  are  constrained  by  the  required  computation  time.  It  is  recommended 
that  the  feasibility  of  the  following  be  investigated  to  decrease  computation  time  : 

•  optimization  of  the  computer  code  at  the  expense  of  program  modularity, 

•  parallel  processing  techniques, 

•  availability  of  a  faster  computing  system  than  the  IBM  3033/438 1  Network. 

The  very  nature  of  the  full-wave  solution  to  the  pulse  propagation  problem  makes  our 
model  computation  intensive.  If  this  method  of  solution  is  to  be  a  useful  tool  in 
investigating  realistic  scenarios,  we  need  to  consider  either  increasing  the  speed  with  which 
the  current  calculations  are  handled,  or  use  another  technique,  such  as  the  Fast-Field- 
Program  method  of  evaluating  the  wave  vector  component  integrals. 
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